diff --git a/plugins/MRCC_Utils/amplitudes.irp.f b/plugins/MRCC_Utils/amplitudes.irp.f index 718d5340..2694aa75 100644 --- a/plugins/MRCC_Utils/amplitudes.irp.f +++ b/plugins/MRCC_Utils/amplitudes.irp.f @@ -137,6 +137,7 @@ END_PROVIDER BEGIN_PROVIDER [ integer, mrcc_AtA_ind, (N_det_ref * n_exc_active) ] &BEGIN_PROVIDER [ double precision, mrcc_AtA_val, (N_states, N_det_ref * n_exc_active) ] +&BEGIN_PROVIDER [ double precision, mrcc_AtS2A_val, (N_states, N_det_ref * n_exc_active) ] &BEGIN_PROVIDER [ integer, mrcc_col_shortcut, (n_exc_active) ] &BEGIN_PROVIDER [ integer, mrcc_N_col, (n_exc_active) ] implicit none @@ -145,11 +146,15 @@ END_PROVIDER END_DOC integer :: AtA_size, i,k integer :: at_roww, at_row, wk, a_coll, a_col, r1, r2, s - double precision, allocatable :: t(:), A_val_mwen(:,:) + double precision, allocatable :: t(:), ts(:), A_val_mwen(:,:), As2_val_mwen(:,:) integer, allocatable :: A_ind_mwen(:) + integer(bit_kind) :: det1(N_int,2), det2(N_int,2) + double precision :: sij + PROVIDE psi_non_ref S_z2_Sz S_z mrcc_AtA_ind(:) = 0 mrcc_AtA_val(:,:) = 0.d0 + mrcc_AtS2A_val(:,:) = 0.d0 mrcc_col_shortcut(:) = 0 mrcc_N_col(:) = 0 AtA_size = 0 @@ -157,9 +162,11 @@ END_PROVIDER !$OMP PARALLEL default(none) shared(k, active_excitation_to_determinants_idx,& !$OMP active_excitation_to_determinants_val, hh_nex) & - !$OMP private(at_row, a_col, t, i, r1, r2, wk, A_ind_mwen, A_val_mwen, a_coll, at_roww)& - !$OMP shared(N_states,mrcc_col_shortcut, mrcc_N_col, AtA_size, mrcc_AtA_val, mrcc_AtA_ind, n_exc_active, active_pp_idx) - allocate(A_val_mwen(N_states,hh_nex), A_ind_mwen(hh_nex), t(N_states)) + !$OMP private(at_row, a_col, t, ts, i, r1, r2, wk, A_ind_mwen, A_val_mwen,& + !$OMP det1,det2,As2_val_mwen, a_coll, at_roww,sij) & + !$OMP shared(N_states,mrcc_col_shortcut, mrcc_N_col, AtA_size, mrcc_AtA_val, mrcc_AtA_ind, & + !$OMP n_exc_active, active_pp_idx,psi_non_ref,N_int,S_z2_Sz, mrcc_AtS2A_val) + allocate(A_val_mwen(N_states,hh_nex), As2_val_mwen(N_states,hh_nex), A_ind_mwen(hh_nex), t(N_states), ts(N_states)) !$OMP DO schedule(dynamic, 100) do at_roww = 1, n_exc_active ! hh_nex @@ -170,6 +177,7 @@ END_PROVIDER do a_coll = 1, n_exc_active a_col = active_pp_idx(a_coll) t(:) = 0d0 + ts(:) = 0d0 r1 = 1 r2 = 1 do while ((active_excitation_to_determinants_idx(r1, at_roww) /= 0).and.(active_excitation_to_determinants_idx(r2, a_coll) /= 0)) @@ -178,8 +186,12 @@ END_PROVIDER else if(active_excitation_to_determinants_idx(r1, at_roww) < active_excitation_to_determinants_idx(r2, a_coll)) then r1 = r1+1 else + det1(:,:) = psi_non_ref(:,:, active_excitation_to_determinants_idx(r1,at_roww)) + det2(:,:) = psi_non_ref(:,:, active_excitation_to_determinants_idx(r2,a_coll)) + call get_s2(det1, det2,N_int,sij) do s=1,N_states t(s) = t(s) - active_excitation_to_determinants_val(s,r1, at_roww) * active_excitation_to_determinants_val(s,r2, a_coll) + ts(s) = ts(s) - active_excitation_to_determinants_val(s,r1, at_roww) * active_excitation_to_determinants_val(s,r2, a_coll) * sij enddo r1 = r1+1 r2 = r2+1 @@ -189,6 +201,7 @@ END_PROVIDER if(a_col == at_row) then do s=1,N_states t(s) = t(s) + 1.d0 + ts(s) = ts(s) + S_z2_Sz enddo end if if(sum(abs(t)) /= 0.d0) then @@ -196,6 +209,7 @@ END_PROVIDER A_ind_mwen(wk) = a_col do s=1,N_states A_val_mwen(s,wk) = t(s) + As2_val_mwen(s,wk) = ts(s) enddo end if end do @@ -212,6 +226,7 @@ END_PROVIDER mrcc_AtA_ind(AtA_size+i) = A_ind_mwen(i) do s=1,N_states mrcc_AtA_val(s,AtA_size+i) = A_val_mwen(s,i) + mrcc_AtS2A_val(s,AtA_size+i) = As2_val_mwen(s,i) enddo enddo AtA_size += wk @@ -219,7 +234,7 @@ END_PROVIDER end if end do !$OMP END DO NOWAIT - deallocate (A_ind_mwen, A_val_mwen, t) + deallocate (A_ind_mwen, A_val_mwen, As2_val_mwen, t, ts) !$OMP END PARALLEL print *, "ATA SIZE", ata_size diff --git a/plugins/MRCC_Utils/mrcc_utils.irp.f b/plugins/MRCC_Utils/mrcc_utils.irp.f index 191866aa..f864dd08 100644 --- a/plugins/MRCC_Utils/mrcc_utils.irp.f +++ b/plugins/MRCC_Utils/mrcc_utils.irp.f @@ -626,8 +626,6 @@ END_PROVIDER double precision :: norm, cx, res integer, allocatable :: lref(:), A_ind_mwen(:) double precision :: phase -! double precision , allocatable :: mrcc_AtA_val(:,:) -! integer, allocatable :: mrcc_AtA_ind(:), col_shortcut(:), , mrcc_N_col(:) double precision, allocatable :: rho_mrcc_init(:,:) @@ -635,13 +633,10 @@ END_PROVIDER print *, "TI", hh_nex, N_det_non_ref - - allocate(rho_mrcc_init(N_det_non_ref, N_states)) - + allocate(x_new(hh_nex)) allocate(x(hh_nex), AtB(hh_nex)) x = 0d0 - allocate(x_new(hh_nex)) do s=1,N_states @@ -712,28 +707,37 @@ END_PROVIDER x_new = x + double precision :: s2(N_states), s2_local, dx double precision :: factor, resold factor = 1.d0 resold = huge(1.d0) do k=0,100000 - !$OMP PARALLEL default(shared) private(cx, i, j, a_col, a_coll) + !$OMP PARALLEL default(shared) private(cx, dx, i, j, a_col, a_coll, s2_local) !$OMP DO do i=1,N_det_non_ref - rho_mrcc(i,s) = rho_mrcc_init(i,s) ! 0d0 + rho_mrcc(i,s) = rho_mrcc_init(i,s) enddo !$OMP END DO + s2(s) = 0.d0 !$OMP DO - do a_coll = 1, n_exc_active !: hh_nex + do a_coll = 1, n_exc_active a_col = active_pp_idx(a_coll) - cx = 0d0 + cx = 0.d0 + dx = 0.d0 do i=mrcc_col_shortcut(a_coll), mrcc_col_shortcut(a_coll) + mrcc_N_col(a_coll) - 1 cx = cx + x(mrcc_AtA_ind(i)) * mrcc_AtA_val(s,i) + dx = dx + x(mrcc_AtA_ind(i)) * mrcc_AtS2A_val(s,i) + s2_local = s2_local + X(a_col)*X(mrcc_AtA_ind(i))*mrcc_AtS2A_val(s,i) end do - x_new(a_col) = AtB(a_col) + cx * factor + x_new(a_col) = AtB(a_col) + (cx+dx) * factor end do !$OMP END DO + + !$OMP CRITICAL + s2(s) = s2(s) + s2_local + !$OMP END CRITICAL !$OMP END PARALLEL @@ -756,14 +760,13 @@ END_PROVIDER resold = res if(mod(k, 100) == 0) then - print *, "res ", k, res, factor + print *, "res ", k, res, s2(s) end if if(res < 1d-9) exit end do - norm = 0.d0 do i=1,N_det_non_ref norm = norm + rho_mrcc(i,s)*rho_mrcc(i,s) diff --git a/src/Determinants/s2.irp.f b/src/Determinants/s2.irp.f index c6bb8390..7e62befb 100644 --- a/src/Determinants/s2.irp.f +++ b/src/Determinants/s2.irp.f @@ -1,36 +1,36 @@ subroutine get_s2(key_i,key_j,Nint,s2) - implicit none - use bitmasks - BEGIN_DOC -! Returns - END_DOC - integer, intent(in) :: Nint - integer(bit_kind), intent(in) :: key_i(Nint,2) - integer(bit_kind), intent(in) :: key_j(Nint,2) - double precision, intent(out) :: s2 - integer :: exc(0:2,2,2) - integer :: degree - double precision :: phase_spsm - integer :: nup, i - - s2 = 0.d0 - !$FORCEINLINE - call get_excitation_degree(key_i,key_j,degree,Nint) - select case (degree) - case(2) - call get_double_excitation(key_j,key_i,exc,phase_spsm,Nint) - if (exc(0,1,1) == 1) then ! Mono alpha + mono-beta - if ( (exc(1,1,1) == exc(1,2,2)).and.(exc(1,1,2) == exc(1,2,1)) ) then - s2 = -phase_spsm - endif - endif - case(0) + implicit none + use bitmasks + BEGIN_DOC + ! Returns + END_DOC + integer, intent(in) :: Nint + integer(bit_kind), intent(in) :: key_i(Nint,2) + integer(bit_kind), intent(in) :: key_j(Nint,2) + double precision, intent(out) :: s2 + integer :: exc(0:2,2,2) + integer :: degree + double precision :: phase_spsm + integer :: nup, i + + s2 = 0.d0 + !$FORCEINLINE + call get_excitation_degree(key_i,key_j,degree,Nint) + select case (degree) + case(2) + call get_double_excitation(key_j,key_i,exc,phase_spsm,Nint) + if (exc(0,1,1) == 1) then ! Mono alpha + mono-beta + if ( (exc(1,1,1) == exc(1,2,2)).and.(exc(1,1,2) == exc(1,2,1)) ) then + s2 = -phase_spsm + endif + endif + case(0) nup = 0 do i=1,Nint nup += popcnt(iand(xor(key_i(i,1),key_i(i,2)),key_i(i,1))) enddo s2 = dble(nup) - end select + end select end BEGIN_PROVIDER [ double precision, S_z ]