10
0
mirror of https://github.com/LCPQ/quantum_package synced 2024-06-01 19:05:25 +02:00

N_states > 1

This commit is contained in:
Yann Garniron 2016-07-06 15:43:21 +02:00
parent e319cb7f1d
commit 6481083bc3
4 changed files with 114 additions and 334 deletions

View File

@ -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/<Psi_0|H|D_m> 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/<Psi_0|H|D_m> 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

View File

@ -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
! <I| <> |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
! <I| /k\ |alpha>
! <I|H|k>
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

View File

@ -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')

View File

@ -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