diff --git a/plugins/MRCC_Utils/mrcc_utils.irp.f b/plugins/MRCC_Utils/mrcc_utils.irp.f index 513e7d09..6c2eb133 100644 --- a/plugins/MRCC_Utils/mrcc_utils.irp.f +++ b/plugins/MRCC_Utils/mrcc_utils.irp.f @@ -1,93 +1,11 @@ use bitmasks BEGIN_PROVIDER [ integer, mrmode ] -&BEGIN_PROVIDER [ logical, old_lambda ] -&BEGIN_PROVIDER [ logical, no_mono_dressing ] - implicit none - CHARACTER(len=255) :: test - CALL get_environment_variable("OLD_LAMBDA", test) - old_lambda = trim(test) /= "" .and. trim(test) /= "0" - CALL get_environment_variable("NO_MONO_DRESSING", test) - no_mono_dressing = trim(test) /= "" .and. trim(test) /= "0" - print *, "old", old_lambda, "mono", no_mono_dressing mrmode = 0 END_PROVIDER - -BEGIN_PROVIDER [ double precision, lambda_mrcc_old, (N_states,psi_det_size) ] -&BEGIN_PROVIDER [ integer, lambda_mrcc_pt2_old, (0:psi_det_size) ] -&BEGIN_PROVIDER [ integer, lambda_mrcc_pt3_old, (0:psi_det_size) ] - implicit none - BEGIN_DOC - cm/ or perturbative 1/Delta_E(m) - END_DOC - integer :: i,k - double precision :: ihpsi_current(N_states) - integer :: i_pert_count - double precision :: hii, lambda_pert - integer :: N_lambda_mrcc_pt2, N_lambda_mrcc_pt3 - double precision, parameter :: x = 2.d0 - double precision :: nurm - i_pert_count = 0 - lambda_mrcc_old = 0.d0 - N_lambda_mrcc_pt2 = 0 - N_lambda_mrcc_pt3 = 0 - lambda_mrcc_pt2_old(0) = 0 - lambda_mrcc_pt3_old(0) = 0 - if(N_states > 1) stop "old lambda N_states == 1" - nurm = 0d0 - do i=1,N_det_ref - nurm += psi_ref_coef(i,1)**2 - end do - - do i=1,N_det_non_ref - call i_h_psi(psi_non_ref(1,1,i), psi_ref, psi_ref_coef, N_int, N_det_ref, & - size(psi_ref_coef,1), N_states,ihpsi_current) - call i_H_j(psi_non_ref(1,1,i),psi_non_ref(1,1,i),N_int,hii) - do k=1,N_states - if (ihpsi_current(k) == 0.d0) then - ihpsi_current(k) = 1.d-32 - endif - lambda_mrcc_old(k,i) = psi_non_ref_coef(i,k)/ihpsi_current(k) - !if ( dabs(psi_non_ref_coef(i,k)*ihpsi_current(k)) < 1.d-5 .or. lambda_mrcc_old(k,i) > 0d0) then - if ( dabs(ihpsi_current(k))*sqrt(psi_non_ref_coef(i,k)**2 / nurm) < 1.d-5 .or. lambda_mrcc_old(k,i) > 0d0) then - i_pert_count += 1 - lambda_mrcc_old(k,i) = 0.d0 - if (lambda_mrcc_pt2_old(N_lambda_mrcc_pt2) /= i) then - N_lambda_mrcc_pt2 += 1 - lambda_mrcc_pt2_old(N_lambda_mrcc_pt2) = i - endif - else - if (lambda_mrcc_pt3_old(N_lambda_mrcc_pt3) /= i) then - N_lambda_mrcc_pt3 += 1 - lambda_mrcc_pt3_old(N_lambda_mrcc_pt3) = i - endif - endif -! lambda_pert = 1.d0 / (psi_ref_energy_diagonalized(k)-hii) -! if((ihpsi_current(k) * lambda_pert) < 0.5d0 * psi_non_ref_coef_restart(i,k) ) then -! lambda_mrcc_old(k,i) = 0.d0 -! endif - - if (lambda_mrcc_old(k,i) > x) then - lambda_mrcc_old(k,i) = x - else if (lambda_mrcc_old(k,i) < -x) then - lambda_mrcc_old(k,i) = -x - endif - enddo - enddo - lambda_mrcc_pt2_old(0) = N_lambda_mrcc_pt2 - lambda_mrcc_pt3_old(0) = N_lambda_mrcc_pt3 - - print*,'N_det_non_ref = ',N_det_non_ref - print*,'Number of ignored determinants = ',i_pert_count - print*,'psi_coef_ref_ratio = ',psi_ref_coef(2,1)/psi_ref_coef(1,1) - print*,'lambda min/max = ',maxval(dabs(lambda_mrcc_old)), minval(dabs(lambda_mrcc_old)) - -END_PROVIDER - - - BEGIN_PROVIDER [ double precision, lambda_mrcc, (N_states,psi_det_size) ] + BEGIN_PROVIDER [ double precision, lambda_mrcc, (N_states, N_det_non_ref) ] &BEGIN_PROVIDER [ integer, lambda_mrcc_pt2, (0:psi_det_size) ] &BEGIN_PROVIDER [ integer, lambda_mrcc_pt3, (0:psi_det_size) ] implicit none @@ -99,49 +17,41 @@ END_PROVIDER integer :: i_pert_count double precision :: hii, lambda_pert integer :: N_lambda_mrcc_pt2, N_lambda_mrcc_pt3 - integer :: histo(200), j - histo = 0 - if(old_lambda) then - lambda_mrcc = lambda_mrcc_old - lambda_mrcc_pt2 = lambda_mrcc_pt2_old - lambda_mrcc_pt3 = lambda_mrcc_pt3_old - else - i_pert_count = 0 - lambda_mrcc = 0.d0 - N_lambda_mrcc_pt2 = 0 - N_lambda_mrcc_pt3 = 0 - lambda_mrcc_pt2(0) = 0 - lambda_mrcc_pt3(0) = 0 + i_pert_count = 0 + lambda_mrcc = 0.d0 + N_lambda_mrcc_pt2 = 0 + N_lambda_mrcc_pt3 = 0 + lambda_mrcc_pt2(0) = 0 + lambda_mrcc_pt3(0) = 0 - do i=1,N_det_non_ref - call i_h_psi(psi_non_ref(1,1,i), psi_ref, psi_ref_coef, N_int, N_det_ref,& - size(psi_ref_coef,1), N_states,ihpsi_current) - call i_H_j(psi_non_ref(1,1,i),psi_non_ref(1,1,i),N_int,hii) - do k=1,N_states - if (ihpsi_current(k) == 0.d0) then - ihpsi_current(k) = 1.d-32 + do i=1,N_det_non_ref + call i_h_psi(psi_non_ref(1,1,i), psi_ref, psi_ref_coef, N_int, N_det_ref,& + size(psi_ref_coef,1), N_states,ihpsi_current) + call i_H_j(psi_non_ref(1,1,i),psi_non_ref(1,1,i),N_int,hii) + do k=1,N_states + if (ihpsi_current(k) == 0.d0) then + ihpsi_current(k) = 1.d-32 + endif + lambda_mrcc(k,i) = min(-1.d-32,psi_non_ref_coef(i,k)/ihpsi_current(k) ) + lambda_pert = 1.d0 / (psi_ref_energy_diagonalized(k)-hii) + if (lambda_pert / lambda_mrcc(k,i) < 0.5d0) then + i_pert_count += 1 + lambda_mrcc(k,i) = 0.d0 + if (lambda_mrcc_pt2(N_lambda_mrcc_pt2) /= i) then + N_lambda_mrcc_pt2 += 1 + lambda_mrcc_pt2(N_lambda_mrcc_pt2) = i endif - lambda_mrcc(k,i) = min(-1.d-32,psi_non_ref_coef(i,k)/ihpsi_current(k) ) - lambda_pert = 1.d0 / (psi_ref_energy_diagonalized(k)-hii) - if (lambda_pert / lambda_mrcc(k,i) < 0.5d0) then - i_pert_count += 1 - lambda_mrcc(k,i) = 0.d0 - if (lambda_mrcc_pt2(N_lambda_mrcc_pt2) /= i) then - N_lambda_mrcc_pt2 += 1 - lambda_mrcc_pt2(N_lambda_mrcc_pt2) = i - endif - else - if (lambda_mrcc_pt3(N_lambda_mrcc_pt3) /= i) then - N_lambda_mrcc_pt3 += 1 - lambda_mrcc_pt3(N_lambda_mrcc_pt3) = i - endif + else + if (lambda_mrcc_pt3(N_lambda_mrcc_pt3) /= i) then + N_lambda_mrcc_pt3 += 1 + lambda_mrcc_pt3(N_lambda_mrcc_pt3) = i endif - enddo + endif enddo - lambda_mrcc_pt2(0) = N_lambda_mrcc_pt2 - lambda_mrcc_pt3(0) = N_lambda_mrcc_pt3 - end if + enddo + lambda_mrcc_pt2(0) = N_lambda_mrcc_pt2 + lambda_mrcc_pt3(0) = N_lambda_mrcc_pt3 print*,'N_det_non_ref = ',N_det_non_ref print*,'psi_coef_ref_ratio = ',psi_ref_coef(2,1)/psi_ref_coef(1,1) print*,'lambda max = ',maxval(dabs(lambda_mrcc)) @@ -149,44 +59,6 @@ END_PROVIDER END_PROVIDER -! BEGIN_PROVIDER [ double precision, lambda_mrcc, (N_states,psi_det_size) ] -! &BEGIN_PROVIDER [ integer, lambda_mrcc_pt2, (0:psi_det_size) ] -! &BEGIN_PROVIDER [ integer, lambda_mrcc_pt3, (0:psi_det_size) ] -! implicit none -! BEGIN_DOC -! ! cm/ or perturbative 1/Delta_E(m) -! END_DOC -! integer :: i,ii,k -! double precision :: ihpsi_current(N_states) -! integer :: i_pert_count -! double precision :: hii, lambda_pert, phase -! integer :: N_lambda_mrcc_pt2, N_lambda_mrcc_pt3, degree -! integer :: exc(N_int, 2) -! histo = 0 -! -! i_pert_count = 0 -! lambda_mrcc = 0.d0 -! N_lambda_mrcc_pt2 = 0 -! N_lambda_mrcc_pt3 = 0 -! lambda_mrcc_pt2(0) = 0 -! lambda_mrcc_pt3(0) = 0 -! -! do ii=1, N_det_ref -! do i=1,N_det_non_ref -! call get_excitation(psi_ref(1,1,II), psi_non_ref(1,1,i), exc, degree, phase, N_int) -! if(degree == -1) cycle -! call i_H_j(psi_non_ref(1,1,ii),psi_non_ref(1,1,i),N_int,hii) -! -! -! lambda_mrcc_pt2(0) = N_lambda_mrcc_pt2 -! lambda_mrcc_pt3(0) = N_lambda_mrcc_pt3 -! -! print*,'N_det_non_ref = ',N_det_non_ref -! print*,'psi_coef_ref_ratio = ',psi_ref_coef(2,1)/psi_ref_coef(1,1) -! print*,'lambda max = ',maxval(dabs(lambda_mrcc)) -! print*,'Number of ignored determinants = ',i_pert_count -! -! END_PROVIDER BEGIN_PROVIDER [ double precision, hij_mrcc, (N_det_non_ref,N_det_ref) ] @@ -362,16 +234,6 @@ logical function is_generable(det1, det2, Nint) return end if if(degree > 2) stop "?22??" - !!!!! -! call dec_exc(exc, h1, h2, p1, p2) -! f = searchExc(toutmoun(1,1), (/h1, h2, p1, p2/), hh_shortcut(hh_shortcut(0)+1)-1) -! !print *, toutmoun(:,1), hh_shortcut(hh_shortcut(0)+1)-1, (/h1, h2, p1, p2/) -! if(f /= -1) then -! is_generable = .true. -! if(.not. excEq(toutmoun(1,f), (/h1, h2, p1, p2/))) stop "????" -! end if -! ! print *, f -! return call decode_exc_int2(exc,degree,h1,p1,h2,p2,s1,s2) @@ -680,10 +542,10 @@ END_PROVIDER END_PROVIDER -BEGIN_PROVIDER [ double precision, dIj, (hh_shortcut(hh_shortcut(0)+1)-1) ] +BEGIN_PROVIDER [ double precision, dIj_unique, (hh_shortcut(hh_shortcut(0)+1)-1, N_states) ] implicit none logical :: ok - integer :: i, j, k, II, pp, hh, ind, wk, nex, a_col, at_row + integer :: i, j, k, s, II, pp, hh, ind, wk, nex, a_col, at_row integer, external :: searchDet, unsortedSearchDet integer(bit_kind) :: myDet(N_int, 2), myMask(N_int, 2) integer :: N, INFO, AtA_size, r1, r2 @@ -691,22 +553,36 @@ BEGIN_PROVIDER [ double precision, dIj, (hh_shortcut(hh_shortcut(0)+1)-1) ] double precision :: t, norm, cx integer, allocatable :: A_ind(:,:), lref(:), AtA_ind(:), A_ind_mwen(:), col_shortcut(:), N_col(:) - if(n_states /= 1) stop "n_states /= 1" + nex = hh_shortcut(hh_shortcut(0)+1)-1 print *, "TI", nex, N_det_non_ref allocate(A_ind(N_det_ref+1, nex), A_val(N_det_ref+1, nex)) - allocate(AtA_ind(N_det_ref * nex), AtA_val(N_det_ref * nex)) !!!!! MAY BE TOO SMALL !!!!!!!! + allocate(AtA_ind(N_det_ref * nex), AtA_val(N_det_ref * nex)) !!!!! MAY BE TOO SMALL ? !!!!!!!! allocate(x(nex), AtB(nex)) allocate(A_val_mwen(nex), A_ind_mwen(nex)) allocate(N_col(nex), col_shortcut(nex), B(N_det_non_ref)) + allocate (x_new(nex)) + + do s = 1, N_states + A_val = 0d0 A_ind = 0 + AtA_ind = 0 + AtA_val = 0d0 + x = 0d0 + AtB = 0d0 + A_val_mwen = 0d0 + A_ind_mwen = 0 + N_col = 0 + col_shortcut = 0 + B = 0d0 + x_new = 0d0 + !$OMP PARALLEL DO schedule(static,10) default(none) shared(psi_non_ref, hh_exists, pp_exists, N_int, A_val, A_ind) & - !$OMP shared(hh_shortcut, psi_ref_coef, N_det_non_ref, psi_non_ref_sorted, psi_non_ref_sorted_idx, psi_ref, N_det_ref) & + !$OMP shared(s, hh_shortcut, psi_ref_coef, N_det_non_ref, psi_non_ref_sorted, psi_non_ref_sorted_idx, psi_ref, N_det_ref) & !$OMP private(lref, pp, II, ok, myMask, myDet, ind, wk) do hh = 1, hh_shortcut(0) - !print *, hh, "/", hh_shortcut(0) do pp = hh_shortcut(hh), hh_shortcut(hh+1)-1 allocate(lref(N_det_non_ref)) lref = 0 @@ -715,12 +591,8 @@ BEGIN_PROVIDER [ double precision, dIj, (hh_shortcut(hh_shortcut(0)+1)-1) ] if(.not. ok) cycle call apply_particle(myMask, pp_exists(1, pp), myDet, ok, N_int) if(.not. ok) cycle - !ind = unsortedSearchDet(psi_non_ref(1,1,1), myDet, N_det_non_ref, N_int) ind = searchDet(psi_non_ref_sorted(1,1,1), myDet(1,1), N_det_non_ref, N_int) if(ind /= -1) then - !iwk = wk+1 - !A_val(wk, pp) = psi_ref_coef(II, 1) - !A_ind(wk, pp) = psi_non_ref_sorted_idx(ind) lref(psi_non_ref_sorted_idx(ind)) = II end if end do @@ -728,7 +600,7 @@ BEGIN_PROVIDER [ double precision, dIj, (hh_shortcut(hh_shortcut(0)+1)-1) ] do i=1, N_det_non_ref if(lref(i) /= 0) then wk += 1 - A_val(wk, pp) = psi_ref_coef(lref(i), 1) + A_val(wk, pp) = psi_ref_coef(lref(i), s) A_ind(wk, pp) = i end if end do @@ -744,13 +616,13 @@ BEGIN_PROVIDER [ double precision, dIj, (hh_shortcut(hh_shortcut(0)+1)-1) ] N_col = 0 !$OMP PARALLEL DO schedule(dynamic, 100) default(none) shared(k, psi_non_ref_coef, A_ind, A_val, x, N_det_ref, nex, N_det_non_ref) & !$OMP private(at_row, a_col, t, i, r1, r2, wk, A_ind_mwen, A_val_mwen) & - !$OMP shared(col_shortcut, N_col, AtB, AtA_size, AtA_val, AtA_ind) + !$OMP shared(col_shortcut, N_col, AtB, AtA_size, AtA_val, AtA_ind, s) do at_row = 1, nex wk = 0 if(mod(at_row, 10000) == 0) print *, "AtA", at_row, "/", nex do i=1,N_det_ref if(A_ind(i, at_row) == 0) exit - AtB(at_row) = AtB(at_row) + psi_non_ref_coef(A_ind(i, at_row), 1) * A_val(i, at_row) + AtB(at_row) = AtB(at_row) + psi_non_ref_coef(A_ind(i, at_row), s) * A_val(i, at_row) end do do a_col = 1, nex t = 0d0 @@ -769,15 +641,11 @@ BEGIN_PROVIDER [ double precision, dIj, (hh_shortcut(hh_shortcut(0)+1)-1) ] end do if(a_col == at_row) then - t = (t + 1d0)! / 2d0 - !print *, a_col, t-1d0 + t = (t + 1d0) end if if(t /= 0d0) then wk += 1 - !AtA_ind(1, wk) = at_row - !AtA_ind(2, wk) = a_col A_ind_mwen(wk) = a_col - !AtA_val(wk) = t A_val_mwen(wk) = t end if end do @@ -796,7 +664,6 @@ BEGIN_PROVIDER [ double precision, dIj, (hh_shortcut(hh_shortcut(0)+1)-1) ] x = AtB if(AtA_size > size(AtA_val)) stop "SIZA" print *, "ATA SIZE", ata_size - allocate (x_new(nex)) integer :: iproc, omp_get_thread_num iproc = omp_get_thread_num() do i=1,nex @@ -821,7 +688,7 @@ BEGIN_PROVIDER [ double precision, dIj, (hh_shortcut(hh_shortcut(0)+1)-1) ] double precision :: norm_cas norm_cas = 0d0 do i = 1, N_det_ref - norm_cas += psi_ref_coef(i,1)**2 + norm_cas += psi_ref_coef(i,s)**2 end do norm = 0d0 @@ -831,23 +698,6 @@ BEGIN_PROVIDER [ double precision, dIj, (hh_shortcut(hh_shortcut(0)+1)-1) ] t = t + X_new(j) * X_new(j) end do - !t = (1d0 - norm_cas) / t - !x_new = x_new * sqrt(t) - !!!!!!!!!!!!!! - !B = 0d0 - !do i=1, nex - ! do j=1, N_det_ref - ! if(A_ind(j, i) == 0) exit - ! B(A_ind(j, i)) += A_val(j, i) * x(i) - ! end do - !end do - !t = 0d0 - !do i=1, size(B) - ! t += B(i)**2 - !end do - !print *, "NORMT", sqrt(t + norm_cas) - !x_new = x_new / sqrt(t + norm_cas) -!!!!!!!!!! t = (1d0 / norm_cas - 1d0) / t x_new = x_new * sqrt(t) @@ -858,7 +708,7 @@ BEGIN_PROVIDER [ double precision, dIj, (hh_shortcut(hh_shortcut(0)+1)-1) ] end do - if(mod(k, 50) == 0) then + if(mod(k, 100) == 0) then print *, "residu ", k, norm, "norm t", sqrt(t) end if @@ -866,77 +716,50 @@ BEGIN_PROVIDER [ double precision, dIj, (hh_shortcut(hh_shortcut(0)+1)-1) ] end do print *, "CONVERGENCE : ", norm + dIj_unique(:size(X), s) = X(:) + -!do k=0,500 -! if(k == 1) print *, X(:10) -! x_new = 0d0 -! A_dense = 0d0 -! !!$OMP PARALLEL DO schedule(dynamic, 10) default(none) shared(k, psi_non_ref_coef, x_new, A_ind, A_val, x, N_det_ref, nex, N_det_non_ref) & -! !!$OMP private(a_col, t, i, cx) & -! !!$OMP firstprivate(A_dense) -! do at_row = 1, nex -! ! ! d DIR$ IVDEP -! cx = 0d0 -! do i=1,N_det_ref -! if(A_ind(i, at_row) == 0) exit -! if(k /= 0) A_dense(A_ind(i, at_row)) = A_val(i, at_row) -! cx = cx + psi_non_ref_coef(A_ind(i, at_row), 1) * A_val(i, at_row) -! !x_new(at_row) = x_new(at_row) + psi_non_ref_coef(A_ind(i, at_row), 1) * A_val(i, at_row) -! end do -! if(k == 0) then -! x_new(at_row) = cx -! cycle -! end if -! do a_col = 1, nex -! t = 0d0 -! do i = 1, N_det_ref -! if(A_ind(i, a_col) == 0) exit -! t = t - A_val(i, a_col) * A_dense(A_ind(i, a_col)) ! -= pcq I-A -! end do -! if(a_col == at_row) t = t + 1d0 -! cx = cx + t * x(a_col) -! !x_new(at_row) = x_new(at_row) + t * x(a_col) -! end do -! x_new(at_row) = cx -! do i=1,N_det_ref -! if(A_ind(i, at_row) == 0) exit -! A_dense(A_ind(i, at_row)) = 0d0 -! end do -! end do -! !!$OMP END PARALLEL DO + end do -! norm = 0d0 -! do j=1, size(X) -! norm += (X_new(j) - X(j))**2 -! X(j) = X_new(j) -! end do -! print *, "residu ", k, norm -! if(norm < 1d-10) exit -!end do -! - dIj(:size(X)) = X(:) - !print *, X print *, "done" END_PROVIDER -double precision function get_dij_index(II, i, Nint) - integer, intent(in) :: II, i, Nint - double precision, external :: get_dij +BEGIN_PROVIDER [ double precision, dij, (N_det_ref, N_det_non_ref, N_states) ] + integer :: s,i,j + print *, "computing amplitudes..." + do s=1, N_states + do i=1, N_det_non_ref + do j=1, N_det_ref + dij(j, i, s) = get_dij_index(j, i, s, N_int) + end do + end do + end do + print *, "done computing amplitudes" +END_PROVIDER - if(dabs(psi_ref_coef(II, 1)) > 1d-1) then - get_dij_index = psi_non_ref_coef(i, 1) / psi_ref_coef(II, 1) + + + +double precision function get_dij_index(II, i, s, Nint) + integer, intent(in) :: II, i, s, Nint + double precision, external :: get_dij + double precision :: HIi + + if(lambda_type == 0) then + get_dij_index = get_dij(psi_ref(1,1,II), psi_non_ref(1,1,i), s, Nint) else - get_dij_index = get_dij(psi_ref(1,1,II), psi_non_ref(1,1,i), Nint) + call i_h_j(psi_ref(1,1,II), psi_non_ref(1,1,i), Nint, HIi) + get_dij_index = HIi * lambda_mrcc(s, i) end if end function -double precision function get_dij(det1, det2, Nint) +double precision function get_dij(det1, det2, s, Nint) use bitmasks implicit none - integer, intent(in) :: Nint + integer, intent(in) :: s, Nint integer(bit_kind) :: det1(Nint, 2), det2(Nint, 2) integer :: degree, f, exc(0:2, 2, 2), t integer*2 :: h1, h2, p1, p2, s1, s2 @@ -976,7 +799,7 @@ double precision function get_dij(det1, det2, Nint) end if if(t /= -1) then - get_dij = dIj(t - 1 + hh_shortcut(f)) + get_dij = dIj_unique(t - 1 + hh_shortcut(f), s) end if end function diff --git a/plugins/mrcepa0/dressing.irp.f b/plugins/mrcepa0/dressing.irp.f index e4b63208..3a91f42e 100644 --- a/plugins/mrcepa0/dressing.irp.f +++ b/plugins/mrcepa0/dressing.irp.f @@ -6,7 +6,7 @@ use bitmasks &BEGIN_PROVIDER [ double precision, delta_ii_mrcc, (N_states, N_det_ref) ] use bitmasks implicit none - integer :: gen, h, p, i_state, n, t, i, h1, h2, p1, p2, s1, s2, iproc + integer :: gen, h, p, n, t, i, h1, h2, p1, p2, s1, s2, iproc integer(bit_kind) :: mask(N_int, 2), omask(N_int, 2) integer(bit_kind),allocatable :: buf(:,:,:) logical :: ok @@ -14,16 +14,16 @@ use bitmasks delta_ij_mrcc = 0d0 delta_ii_mrcc = 0d0 - i_state = 1 + print *, "Dij", dij(1,1,1) provide hh_shortcut psi_det_size! lambda_mrcc !$OMP PARALLEL DO default(none) schedule(dynamic) & !$OMP shared(psi_det_generators, N_det_generators, hh_exists, pp_exists, N_int, hh_shortcut) & - !$OMP shared(N_states, N_det_non_ref, N_det_ref, delta_ii_mrcc, delta_ij_mrcc) & + !$OMP shared(N_det_non_ref, N_det_ref, delta_ii_mrcc, delta_ij_mrcc) & !$OMP private(h, n, mask, omask, buf, ok, iproc) do gen= 1, N_det_generators allocate(buf(N_int, 2, N_det_non_ref)) iproc = omp_get_thread_num() + 1 - print *, gen, "/", N_det_generators + if(mod(gen, 10) == 0) print *, "mrcc ", gen, "/", N_det_generators do h=1, hh_shortcut(0) call apply_hole(psi_det_generators(1,1,gen), hh_exists(1, h), mask, ok, N_int) if(.not. ok) cycle @@ -36,7 +36,9 @@ use bitmasks if(n > N_det_non_ref) stop "MRCC..." end do n = n - 1 + if(n /= 0) call mrcc_part_dress(delta_ij_mrcc, delta_ii_mrcc,gen,n,buf,N_int,omask) + end do deallocate(buf) end do @@ -86,7 +88,8 @@ subroutine mrcc_part_dress(delta_ij_, delta_ii_,i_generator,n_selected,det_buffe integer, allocatable :: idx_microlist(:), N_microlist(:), ptr_microlist(:), idx_microlist_zero(:) integer :: mobiles(2), smallerlist logical, external :: detEq, is_generable - double precision, external :: get_dij, get_dij_index + !double precision, external :: get_dij, get_dij_index + leng = max(N_det_generators, N_det_non_ref) allocate(miniList(Nint, 2, leng), tq(Nint,2,n_selected), idx_minilist(leng), hij_cache(N_det_non_ref)) @@ -171,7 +174,6 @@ subroutine mrcc_part_dress(delta_ij_, delta_ii_,i_generator,n_selected,det_buffe idx_alpha(j) = idx_microlist_zero(idx_alpha(j)) end do - else call get_excitation_degree_vector(miniList,tq(1,1,i_alpha),degree_alpha,Nint,N_minilist,idx_alpha) do j=1,idx_alpha(0) @@ -184,7 +186,6 @@ subroutine mrcc_part_dress(delta_ij_, delta_ii_,i_generator,n_selected,det_buffe 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 @@ -199,7 +200,6 @@ subroutine mrcc_part_dress(delta_ij_, delta_ii_,i_generator,n_selected,det_buffe ! |alpha> do k_sd=1,idx_alpha(0) - ! Loop if lambda == 0 logical :: loop ! loop = .True. @@ -220,18 +220,16 @@ subroutine mrcc_part_dress(delta_ij_, delta_ii_,i_generator,n_selected,det_buffe ! ! - hIk = hij_mrcc(idx_alpha(k_sd),i_I) + !hIk = hij_mrcc(idx_alpha(k_sd),i_I) ! call i_h_j(psi_ref(1,1,i_I),psi_non_ref(1,1,idx_alpha(k_sd)),Nint,hIk) - do i_state=1,N_states - dIK(i_state) = get_dij_index(i_I, idx_alpha(k_sd), Nint) + dIK(i_state) = dij(i_I, idx_alpha(k_sd), i_state) !dIk(i_state) = get_dij(psi_ref(1,1,i_I), psi_non_ref(1,1,idx_alpha(k_sd)), N_int) !!hIk * lambda_mrcc(i_state,idx_alpha(k_sd)) !dIk(i_state) = psi_non_ref_coef(idx_alpha(k_sd), i_state) / psi_ref_coef(i_I, i_state) enddo - ! |l> = Exc(k -> alpha) |I> call get_excitation(psi_non_ref(1,1,idx_alpha(k_sd)),tq(1,1,i_alpha),exc,degree,phase,Nint) call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2) @@ -239,7 +237,6 @@ subroutine mrcc_part_dress(delta_ij_, delta_ii_,i_generator,n_selected,det_buffe tmp_det(k,1) = psi_ref(k,1,i_I) tmp_det(k,2) = psi_ref(k,2,i_I) enddo - logical :: ok call apply_excitation(psi_ref(1,1,i_I), exc, tmp_det, ok, Nint) if(.not. ok) cycle @@ -249,7 +246,6 @@ subroutine mrcc_part_dress(delta_ij_, delta_ii_,i_generator,n_selected,det_buffe 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 @@ -266,7 +262,7 @@ subroutine mrcc_part_dress(delta_ij_, delta_ii_,i_generator,n_selected,det_buffe hIl = hij_mrcc(idx_alpha(l_sd),i_I) ! call i_h_j(psi_ref(1,1,i_I),psi_non_ref(1,1,idx_alpha(l_sd)),Nint,hIl) do i_state=1,N_states - dka(i_state) = get_dij_index(i_I, idx_alpha(l_sd), N_int) * phase * phase2 + dka(i_state) = dij(i_I, idx_alpha(l_sd), i_state) * phase * phase2 !dka(i_state) = get_dij(psi_ref(1,1,i_I), psi_non_ref(1,1,idx_alpha(l_sd)), N_int) * phase * phase2 !hIl * lambda_mrcc(i_state,idx_alpha(l_sd)) * phase * phase2 !dka(i_state) = psi_non_ref_coef(idx_alpha(l_sd), i_state) / psi_ref_coef(i_I, i_state) * phase * phase2 enddo @@ -279,7 +275,7 @@ subroutine mrcc_part_dress(delta_ij_, delta_ii_,i_generator,n_selected,det_buffe dIa(i_state) = dIa(i_state) + dIk(i_state) * dka(i_state) enddo enddo - + do i_state=1,N_states ci_inv(i_state) = psi_ref_coef_inv(i_I,i_state) enddo @@ -292,7 +288,6 @@ subroutine mrcc_part_dress(delta_ij_, delta_ii_,i_generator,n_selected,det_buffe enddo enddo call omp_set_lock( psi_ref_lock(i_I) ) - do i_state=1,N_states if(dabs(psi_ref_coef(i_I,i_state)).ge.5.d-5)then do l_sd=1,idx_alpha(0) @@ -546,12 +541,12 @@ END_PROVIDER implicit none integer :: i,j,k double precision :: Hjk, Hki, Hij - double precision, external :: get_dij + !double precision, external :: get_dij integer i_state, degree provide lambda_mrcc dIj do i_state = 1, N_states - !$OMP PARALLEL DO default(none) schedule(dynamic) private(j,k,Hjk,Hki,degree) shared(no_mono_dressing,lambda_mrcc,i_state, N_det_non_ref,psi_ref, psi_non_ref,N_int,delta_cas,N_det_ref) + !$OMP PARALLEL DO default(none) schedule(dynamic) private(j,k,Hjk,Hki,degree) shared(lambda_mrcc,i_state, N_det_non_ref,psi_ref, psi_non_ref,N_int,delta_cas,N_det_ref,dij) do i=1,N_det_ref do j=1,i call get_excitation_degree(psi_ref(1,1,i), psi_ref(1,1,j), degree, N_int) @@ -561,7 +556,7 @@ END_PROVIDER 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 * get_dij(psi_ref(1,1,i), psi_non_ref(1,1,k), N_int) ! * Hki * lambda_mrcc(i_state, k) + delta_cas(i,j,i_state) += Hjk * dij(i, k, i_state) ! * Hki * lambda_mrcc(i_state, k) !print *, Hjk * get_dij(psi_ref(1,1,i), psi_non_ref(1,1,k), N_int), Hki * get_dij(psi_ref(1,1,j), psi_non_ref(1,1,k), N_int) end do delta_cas(j,i,i_state) = delta_cas(i,j,i_state) @@ -659,7 +654,7 @@ end function integer, allocatable :: idx_sorted_bit(:) integer, external :: get_index_in_psi_det_sorted_bit, searchDet logical, external :: is_in_wavefunction, detEq - double precision, external :: get_dij + !double precision, external :: get_dij integer :: II, blok integer*8, save :: notf = 0 @@ -675,7 +670,7 @@ end function enddo ! To provide everything - contrib = get_dij(psi_ref(1,1,1), psi_non_ref(1,1,1), N_int) + contrib = dij(1, 1, 1) do i_state = 1, N_states delta_mrcepa0_ii(:,:) = 0d0 @@ -685,7 +680,7 @@ end function !$OMP private(m,i,II,J,k,degree,myActive,made_hole,made_particle,hjk,contrib) & !$OMP shared(active_sorb, psi_non_ref, psi_non_ref_coef, psi_ref, psi_ref_coef, cepa0_shortcut, det_cepa0_active) & !$OMP shared(N_det_ref, N_det_non_ref,N_int,det_cepa0_idx,lambda_mrcc,det_ref_active, delta_cas) & - !$OMP shared(notf,i_state, sortRef, sortRefIdx) + !$OMP shared(notf,i_state, sortRef, sortRefIdx, dij) do blok=1,cepa0_shortcut(0) do i=cepa0_shortcut(blok), cepa0_shortcut(blok+1)-1 do II=1,N_det_ref @@ -727,7 +722,7 @@ end function call i_h_j(psi_non_ref(1,1,det_cepa0_idx(k)),psi_ref(1,1,J),N_int,HJk) !contrib = delta_cas(II, J, i_state) * HJk * lambda_mrcc(i_state, det_cepa0_idx(k)) - contrib = delta_cas(II, J, i_state) * get_dij(psi_ref(1,1,J), psi_non_ref(1,1,det_cepa0_idx(k)), N_int) + contrib = delta_cas(II, J, i_state) * dij(J, det_cepa0_idx(k), i_state) !$OMP ATOMIC delta_mrcepa0_ij(J, det_cepa0_idx(i), i_state) += contrib diff --git a/plugins/mrcepa0/dressing_slave.irp.f b/plugins/mrcepa0/dressing_slave.irp.f index 7d64aa5e..7c9d7fe0 100644 --- a/plugins/mrcepa0/dressing_slave.irp.f +++ b/plugins/mrcepa0/dressing_slave.irp.f @@ -55,7 +55,7 @@ subroutine mrsc2_dressing_slave(thread,iproc) logical, external :: is_in_wavefunction, isInCassd, detEq integer,allocatable :: komon(:) logical :: komoned - double precision, external :: get_dij + !double precision, external :: get_dij zmq_to_qp_run_socket = new_zmq_to_qp_run_socket() zmq_socket_push = new_zmq_push_socket(thread) @@ -144,7 +144,7 @@ subroutine mrsc2_dressing_slave(thread,iproc) ! if(I_i == J) phase_Ii = phase_Ji do i_state = 1,N_states - dkI = h_(J,i) * get_dij(psi_ref(1,1,i_I), psi_non_ref(1,1,i), N_int) + dkI = h_(J,i) * dij(i_I, i, i_state)!get_dij(psi_ref(1,1,i_I), psi_non_ref(1,1,i), N_int) !dkI = h_(J,i) * h_(i_I,i) * lambda_mrcc(i_state, i) dleat(i_state, kn, 1) = dkI dleat(i_state, kn, 2) = dkI @@ -174,7 +174,7 @@ subroutine mrsc2_dressing_slave(thread,iproc) !contrib = h_(i_I,k) * lambda_mrcc(i_state, k) * dleat(i_state, m, 2)! * phase_al - contrib = get_dij(psi_ref(1,1,i_I), psi_non_ref(1,1,k), N_int) * dleat(i_state, m, 2) + contrib = dij(i_I, k, i_state) * dleat(i_state, m, 2) delta(i_state,ll,1) += contrib if(dabs(psi_ref_coef(i_I,i_state)).ge.5.d-5) then delta(i_state,0,1) -= contrib * ci_inv(i_state) * psi_non_ref_coef(l,i_state) @@ -182,7 +182,7 @@ subroutine mrsc2_dressing_slave(thread,iproc) if(I_i == J) cycle !contrib = h_(J,l) * lambda_mrcc(i_state, l) * dleat(i_state, m, 1)! * phase_al - contrib = get_dij(psi_ref(1,1,J), psi_non_ref(1,1,l), N_int) * dleat(i_state, m, 1) + contrib = dij(J, l, i_state) * dleat(i_state, m, 1) delta(i_state,kk,2) += contrib if(dabs(psi_ref_coef(J,i_state)).ge.5.d-5) then delta(i_state,0,2) -= contrib * cj_inv(i_state) * psi_non_ref_coef(k,i_state) @@ -483,9 +483,6 @@ end integer :: KKsize = 1000000 - ! -459.6346665282306 - ! -459.6346665282306 - call new_parallel_job(zmq_to_qp_run_socket,'mrsc2') diff --git a/plugins/mrcepa0/mrcepa0_general.irp.f b/plugins/mrcepa0/mrcepa0_general.irp.f index 0ef4c92b..53a0822d 100644 --- a/plugins/mrcepa0/mrcepa0_general.irp.f +++ b/plugins/mrcepa0/mrcepa0_general.irp.f @@ -16,10 +16,11 @@ subroutine run(N_st,energy) double precision :: thresh_mrcc + thresh_mrcc = 1d-7 n_it_mrcc_max = 10 - if(no_mono_dressing) then + if(n_it_mrcc_max == 1) then do j=1,N_states_diag do i=1,N_det psi_coef(i,j) = CI_eigenvectors_dressed(i,j) @@ -73,44 +74,8 @@ subroutine run_pt2(N_st,energy) print*,'Last iteration only to compute the PT2' threshold_selectors = 1.d0 threshold_generators = 0.999d0 - -! N_det_generators = lambda_mrcc_pt2(0) -! do i=1,N_det_generators -! j = lambda_mrcc_pt2(i) -! do k=1,N_int -! psi_det_generators(k,1,i) = psi_non_ref(k,1,j) -! psi_det_generators(k,2,i) = psi_non_ref(k,2,j) -! enddo -! do k=1,N_st -! psi_coef_generators(i,k) = psi_non_ref_coef(j,k) -! enddo -! enddo -! SOFT_TOUCH N_det_generators psi_det_generators psi_coef_generators ci_eigenvectors_dressed ci_eigenvectors_s2_dressed ci_electronic_energy_dressed - - -! -! N_det_generators = lambda_mrcc_pt2(0) + N_det_cas -! do i=1,N_det_cas -! do k=1,N_int -! psi_det_generators(k,1,i) = psi_ref(k,1,i) -! psi_det_generators(k,2,i) = psi_ref(k,2,i) -! enddo -! do k=1,N_st -! psi_coef_generators(i,k) = psi_ref_coef(i,k) -! enddo -! enddo -! do i=N_det_cas+1,N_det_generators -! j = lambda_mrcc_pt2(i - N_det_cas) -! do k=1,N_int -! psi_det_generators(k,1,i) = psi_non_ref(k,1,j) -! psi_det_generators(k,2,i) = psi_non_ref(k,2,j) -! enddo -! do k=1,N_st -! psi_coef_generators(i,k) = psi_non_ref_coef(j,k) -! enddo -! enddo -! SOFT_TOUCH N_det_generators psi_det_generators psi_coef_generators ci_eigenvectors_dressed ci_eigenvectors_s2_dressed ci_electronic_energy_dressed - + + N_det_generators = lambda_mrcc_pt3(0) + N_det_ref N_det_selectors = lambda_mrcc_pt3(0) + N_det_ref