LCPQ/quantum_package
10
0
Fork 0
mirror of https://github.com/LCPQ/quantum_package synced 2024-12-23 12:56:14 +01:00
Code Issues Releases Wiki Activity

Refactoring MRCC

Browse Source
This commit is contained in:
Anthony Scemama 2016-08-02 13:45:17 +02:00
parent eb15a392be
commit 33f5f44fe5
1 changed files with 200 additions and 192 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

108
plugins/MRCC_Utils/mrcc_utils.irp.f
View File

@ -225,6 +225,7 @@ logical function is_generable(det1, det2, Nint)
integer, external :: searchExc integer, external :: searchExc
logical, external :: excEq logical, external :: excEq
double precision :: phase double precision :: phase
integer*2 :: tmp_array(4)
is_generable = .false. is_generable = .false.
call get_excitation(det1, det2, exc, degree, phase, Nint) call get_excitation(det1, det2, exc, degree, phase, Nint)
@ -247,19 +248,20 @@ logical function is_generable(det1, det2, Nint)
end if end if
if(h1 + (s1-1)*mo_tot_num < h2 + (s2-1)*mo_tot_num) then if(h1 + (s1-1)*mo_tot_num < h2 + (s2-1)*mo_tot_num) then
f = searchExc(hh_exists(1,1), (/s1, h1, s2, h2/), hh_shortcut(0)) tmp_array = (/s1, h1, s2, h2/)
else else
f = searchExc(hh_exists(1,1), (/s2, h2, s1, h1/), hh_shortcut(0)) tmp_array = (/s2, h2, s1, h1/)
end if end if
if(f == -1) return f = searchExc(hh_exists(1,1), tmp_array, hh_shortcut(0))
if(p1 + (s1-1)*mo_tot_num < p2 + (s2-1)*mo_tot_num) then if(p1 + (s1-1)*mo_tot_num < p2 + (s2-1)*mo_tot_num) then
f = searchExc(pp_exists(1,hh_shortcut(f)), (/s1, p1, s2, p2/), hh_shortcut(f+1)-hh_shortcut(f)) tmp_array = (/s1, p1, s2, p2/)
else else
f = searchExc(pp_exists(1,hh_shortcut(f)), (/s2, p2, s1, p1/), hh_shortcut(f+1)-hh_shortcut(f)) tmp_array = (/s2, p2, s1, p1/)
end if end if
f = searchExc(pp_exists(1,hh_shortcut(f)), tmp_array, hh_shortcut(f+1)-hh_shortcut(f))
if(f /= -1) is_generable = .true. is_generable = (f /= -1)
end function end function
@ -543,14 +545,15 @@ END_PROVIDER
BEGIN_PROVIDER [ double precision, dIj_unique, (hh_shortcut(hh_shortcut(0)+1)-1, N_states) ] BEGIN_PROVIDER [ double precision, dIj_unique, (hh_shortcut(hh_shortcut(0)+1)-1, N_states) ]
&BEGIN_PROVIDER [ double precision, rho_mrcc, (N_det_non_ref, N_states) ]
implicit none implicit none
logical :: ok logical :: ok
integer :: i, j, k, s, 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, external :: searchDet, unsortedSearchDet
integer(bit_kind) :: myDet(N_int, 2), myMask(N_int, 2) integer(bit_kind) :: myDet(N_int, 2), myMask(N_int, 2)
integer :: N, INFO, AtA_size, r1, r2 integer :: N, INFO, AtA_size, r1, r2
double precision , allocatable:: B(:), AtB(:), AtA_val(:), A_val(:,:), x(:), x_new(:), A_val_mwen(:) double precision , allocatable :: AtB(:), AtA_val(:), A_val(:,:), x(:), x_new(:), A_val_mwen(:)
double precision :: t, norm, cx double precision :: t, norm, cx, res
integer, allocatable :: A_ind(:,:), lref(:), AtA_ind(:), A_ind_mwen(:), col_shortcut(:), N_col(:) integer, allocatable :: A_ind(:,:), lref(:), AtA_ind(:), A_ind_mwen(:), col_shortcut(:), N_col(:)
@ -561,7 +564,7 @@ BEGIN_PROVIDER [ double precision, dIj_unique, (hh_shortcut(hh_shortcut(0)+1)-1,
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(x(nex), AtB(nex))
allocate(A_val_mwen(nex), A_ind_mwen(nex)) allocate(A_val_mwen(nex), A_ind_mwen(nex))
allocate(N_col(nex), col_shortcut(nex), B(N_det_non_ref)) allocate(N_col(nex), col_shortcut(nex))
allocate (x_new(nex)) allocate (x_new(nex))
do s = 1, N_states do s = 1, N_states
@ -571,13 +574,10 @@ BEGIN_PROVIDER [ double precision, dIj_unique, (hh_shortcut(hh_shortcut(0)+1)-1,
AtA_ind = 0 AtA_ind = 0
AtA_val = 0d0 AtA_val = 0d0
x = 0d0 x = 0d0
AtB = 0d0
A_val_mwen = 0d0 A_val_mwen = 0d0
A_ind_mwen = 0 A_ind_mwen = 0
N_col = 0 N_col = 0
col_shortcut = 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 PARALLEL DO schedule(static,10) default(none) shared(psi_non_ref, hh_exists, pp_exists, N_int, A_val, A_ind)&
!$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 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)&
@ -661,66 +661,70 @@ BEGIN_PROVIDER [ double precision, dIj_unique, (hh_shortcut(hh_shortcut(0)+1)-1,
end if end if
end do end do
x = AtB
if(AtA_size > size(AtA_val)) stop "SIZA" if(AtA_size > size(AtA_val)) stop "SIZA"
print *, "ATA SIZE", ata_size print *, "ATA SIZE", ata_size
integer :: iproc, omp_get_thread_num integer :: iproc, omp_get_thread_num
iproc = omp_get_thread_num() iproc = omp_get_thread_num()
do i=1,nex do i=1,nex
x_new(i) = 0.D0 x(i) = AtB(i)
enddo enddo
do k=0,100000 do k=0,100000
!$OMP PARALLEL DO default(shared) !$OMP PARALLEL default(shared) private(cx, i, j, a_col)
do i=1,nex
x_new(i) = AtB(i)
enddo
!$OMP PARALLEL DO default(shared) private(cx, i) !$OMP DO
do i=1,N_det_non_ref
rho_mrcc(i,s) = 0.d0
enddo
!$OMP END DO
!$OMP DO
do a_col = 1, nex do a_col = 1, nex
cx = 0d0 cx = 0d0
do i=col_shortcut(a_col), col_shortcut(a_col) + N_col(a_col) - 1 do i=col_shortcut(a_col), col_shortcut(a_col) + N_col(a_col) - 1
cx += x(AtA_ind(i)) * AtA_val(i) cx = cx + x(AtA_ind(i)) * AtA_val(i)
end do end do
x_new(a_col) += cx x_new(a_col) = AtB(a_col) + cx
end do end do
!$OMP END PARALLEL DO !$OMP END DO
double precision :: norm_cas
norm_cas = 0d0 !$OMP END PARALLEL
do i = 1, N_det_ref
norm_cas += psi_ref_coef(i,s)**2 res = 0.d0
do a_col=1,nex
do j=1,N_det_non_ref
i = A_ind(j,a_col)
if (i==0) exit
rho_mrcc(i,s) = rho_mrcc(i,s) + A_val(j,a_col) * X_new(a_col)
enddo enddo
res = res + (X_new(a_col) - X(a_col))**2
norm = 0d0 X(a_col) = X_new(a_col)
t = 0d0
do j=1, size(X)
t = t + X_new(j) * X_new(j)
end do end do
t = (1d0 - norm_cas ) / t
x_new = x_new * sqrt(t)
do j=1, size(X)
norm += (X_new(j) - X(j))**2
x(j) = x_new(j)
end do
if(mod(k, 100) == 0) then if(mod(k, 100) == 0) then
print *, "residu ", k, norm, "norm t", sqrt(t) print *, "residu ", k, res
end if end if
if(norm < 1d-16) exit if(res < 1d-16) exit
end do end do
print *, "CONVERGENCE : ", norm
norm = 0.d0
do i=1,N_det_non_ref
norm = norm + rho_mrcc(i,s)*rho_mrcc(i,s)
enddo
do i=1,N_det_ref
norm = norm + psi_ref_coef(i,s)*psi_ref_coef(i,s)
enddo
print *, k, "res : ", res, "norm : ", sqrt(norm)
dIj_unique(:size(X), s) = X(:) dIj_unique(:size(X), s) = X(:)
do i=1,N_det_non_ref
rho_mrcc(i,s) = psi_non_ref_coef(i,s) / rho_mrcc(i,s)
enddo enddo
end do
print *, "done" print *, "done"
END_PROVIDER END_PROVIDER
@ -767,6 +771,7 @@ double precision function get_dij(det1, det2, s, Nint)
integer, external :: searchExc integer, external :: searchExc
logical, external :: excEq logical, external :: excEq
double precision :: phase double precision :: phase
integer*2 :: tmp_array(4)
get_dij = 0d0 get_dij = 0d0
call get_excitation(det1, det2, exc, degree, phase, Nint) call get_excitation(det1, det2, exc, degree, phase, Nint)
@ -787,17 +792,20 @@ double precision function get_dij(det1, det2, s, Nint)
end if end if
if(h1 + (s1-1)*mo_tot_num < h2 + (s2-1)*mo_tot_num) then if(h1 + (s1-1)*mo_tot_num < h2 + (s2-1)*mo_tot_num) then
f = searchExc(hh_exists(1,1), (/s1, h1, s2, h2/), hh_shortcut(0)) tmp_array = (/s1, h1, s2, h2/)
else else
f = searchExc(hh_exists(1,1), (/s2, h2, s1, h1/), hh_shortcut(0)) tmp_array = (/s2, h2, s1, h1/)
end if end if
f = searchExc(hh_exists(1,1), tmp_array, hh_shortcut(0))
if(f == -1) return if(f == -1) return
if(p1 + (s1-1)*mo_tot_num < p2 + (s2-1)*mo_tot_num) then if(p1 + (s1-1)*mo_tot_num < p2 + (s2-1)*mo_tot_num) then
t = searchExc(pp_exists(1,hh_shortcut(f)), (/s1, p1, s2, p2/), hh_shortcut(f+1)-hh_shortcut(f)) tmp_array = (/s1, p1, s2, p2/)
else else
t = searchExc(pp_exists(1,hh_shortcut(f)), (/s2, p2, s1, p1/), hh_shortcut(f+1)-hh_shortcut(f)) tmp_array = (/s2, p2, s1, p1/)
end if end if
t = searchExc(pp_exists(1,hh_shortcut(f)), tmp_array, hh_shortcut(f+1)-hh_shortcut(f))
if(t /= -1) then if(t /= -1) then
get_dij = dIj_unique(t - 1 + hh_shortcut(f), s) get_dij = dIj_unique(t - 1 + hh_shortcut(f), s)