From 8a5a671e314199653433c9f01ea7b6074f6ff978 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Fri, 3 Mar 2017 10:41:55 +0100 Subject: [PATCH] Reduced error in MRCC --- plugins/MRCC_Utils/davidson.irp.f | 67 +++++++++-------------------- plugins/MRCC_Utils/mrcc_utils.irp.f | 42 ++++++++++++------ plugins/mrcepa0/dressing.irp.f | 18 +++----- plugins/mrcepa0/mrcc.irp.f | 2 +- 4 files changed, 56 insertions(+), 73 deletions(-) diff --git a/plugins/MRCC_Utils/davidson.irp.f b/plugins/MRCC_Utils/davidson.irp.f index 6bdadb24..c884b3c2 100644 --- a/plugins/MRCC_Utils/davidson.irp.f +++ b/plugins/MRCC_Utils/davidson.irp.f @@ -35,21 +35,20 @@ subroutine davidson_diag_mrcc(dets_in,u_in,energies,dim_in,sze,N_st,N_st_diag,Ni PROVIDE mo_bielec_integrals_in_map allocate(H_jj(sze)) + H_jj(1) = diag_h_mat_elem(dets_in(1,1,1),Nint) !$OMP PARALLEL DEFAULT(NONE) & !$OMP SHARED(sze,H_jj,N_det_ref,dets_in,Nint,istate,delta_ii,idx_ref) & !$OMP PRIVATE(i) - !$OMP DO SCHEDULE(guided) - do i=1,sze + !$OMP DO + do i=2,sze H_jj(i) = diag_h_mat_elem(dets_in(1,1,i),Nint) enddo !$OMP END DO - !$OMP DO SCHEDULE(guided) - do i=1,N_det_ref - H_jj(idx_ref(i)) += delta_ii(istate,i) - enddo - !$OMP END DO !$OMP END PARALLEL + do i=1,N_det_ref + H_jj(idx_ref(i)) += delta_ii(istate,i) + enddo call davidson_diag_hjj_mrcc(dets_in,u_in,H_jj,energies,dim_in,sze,N_st,N_st_diag,Nint,iunit,istate) deallocate (H_jj) end @@ -224,17 +223,6 @@ subroutine davidson_diag_hjj_mrcc(dets_in,u_in,H_jj,energies,dim_in,sze,N_st,N_s W(i,k,iter+1) = 0.d0 enddo enddo -! do k=1,N_st_diag -! do iter2=1,iter -! do l=1,N_st_diag -! do i=1,sze -! U(i,k,iter+1) = U(i,k,iter+1) + U(i,l,iter2)*y(l,iter2,k,1) -! W(i,k,iter+1) = W(i,k,iter+1) + W(i,l,iter2)*y(l,iter2,k,1) -! enddo -! enddo -! enddo -! enddo -! ! call dgemm('N','N', sze, N_st_diag, N_st_diag*iter, & 1.d0, U, size(U,1), y, size(y,1)*size(y,2), 0.d0, U(1,1,iter+1), size(U,1)) @@ -276,27 +264,11 @@ subroutine davidson_diag_hjj_mrcc(dets_in,u_in,H_jj,energies,dim_in,sze,N_st,N_s do k=1,N_st_diag -! do iter2=1,iter -! do l=1,N_st_diag -! c(1) = u_dot_v(U(1,k,iter+1),U(1,l,iter2),sze) -! do i=1,sze -! U(i,k,iter+1) = U(i,k,iter+1) - c(1) * U(i,l,iter2) -! enddo -! enddo -! enddo -! call dgemv('T',sze,N_st_diag*iter,1.d0,U,size(U,1), & U(1,k,iter+1),1,0.d0,c,1) call dgemv('N',sze,N_st_diag*iter,-1.d0,U,size(U,1), & c,1,1.d0,U(1,k,iter+1),1) -! -! do l=1,k-1 -! c(1) = u_dot_v(U(1,k,iter+1),U(1,l,iter+1),sze) -! do i=1,sze -! U(i,k,iter+1) = U(i,k,iter+1) - c(1) * U(i,l,iter+1) -! enddo -! enddo -! + call dgemv('T',sze,k-1,1.d0,U(1,1,iter+1),size(U,1), & U(1,k,iter+1),1,0.d0,c,1) call dgemv('N',sze,k-1,-1.d0,U(1,1,iter+1),size(U,1), & @@ -429,7 +401,7 @@ subroutine H_u_0_mrcc_nstates(v_0,u_0,H_jj,n,keys_tmp,Nint,istate_in,N_st,sze_8) allocate(vt(sze_8,N_st)) Vt = 0.d0 - !$OMP DO SCHEDULE(dynamic) + !$OMP DO SCHEDULE(static,1) do sh=1,shortcut(0,1) do sh2=sh,shortcut(0,1) exa = 0 @@ -468,9 +440,9 @@ subroutine H_u_0_mrcc_nstates(v_0,u_0,H_jj,n,keys_tmp,Nint,istate_in,N_st,sze_8) enddo enddo enddo - !$OMP END DO NOWAIT + !$OMP END DO - !$OMP DO SCHEDULE(dynamic) + !$OMP DO SCHEDULE(static,1) do sh=1,shortcut(0,2) do i=shortcut(sh,2),shortcut(sh+1,2)-1 org_i = sort_idx(i,2) @@ -490,7 +462,7 @@ subroutine H_u_0_mrcc_nstates(v_0,u_0,H_jj,n,keys_tmp,Nint,istate_in,N_st,sze_8) end do end do enddo - !$OMP END DO NOWAIT + !$OMP END DO !$OMP DO do ii=1,n_det_ref @@ -559,25 +531,26 @@ subroutine davidson_diag_mrcc_hs2(dets_in,u_in,dim_in,energies,sze,N_st,N_st_dia ASSERT (sze > 0) ASSERT (Nint > 0) ASSERT (Nint == N_int) - PROVIDE mo_bielec_integrals_in_map + PROVIDE mo_bielec_integrals_in_map allocate(H_jj(sze), S2_jj(sze)) + H_jj(1) = diag_h_mat_elem(dets_in(1,1,1),Nint) + call get_s2(dets_in(1,1,1),dets_in(1,1,1),Nint,S2_jj(1)) !$OMP PARALLEL DEFAULT(NONE) & !$OMP SHARED(sze,H_jj,S2_jj, dets_in,Nint,N_det_ref,delta_ii, & !$OMP idx_ref, istate) & !$OMP PRIVATE(i) - !$OMP DO SCHEDULE(guided) - do i=1,sze + !$OMP DO + do i=2,sze H_jj(i) = diag_h_mat_elem(dets_in(1,1,i),Nint) call get_s2(dets_in(1,1,i),dets_in(1,1,i),Nint,S2_jj(i)) enddo !$OMP END DO - !$OMP DO SCHEDULE(guided) + !$OMP END PARALLEL + do i=1,N_det_ref H_jj(idx_ref(i)) += delta_ii(istate,i) enddo - !$OMP END DO - !$OMP END PARALLEL call davidson_diag_hjj_sjj_mrcc(dets_in,u_in,H_jj,S2_jj,energies,dim_in,sze,N_st,N_st_diag,Nint,iunit,istate) deallocate (H_jj,S2_jj) @@ -1051,7 +1024,7 @@ subroutine H_S2_u_0_mrcc_nstates(v_0,s_0,u_0,H_jj,S2_jj,n,keys_tmp,Nint,istate_i Vt = 0.d0 St = 0.d0 - !$OMP DO SCHEDULE(guided) + !$OMP DO SCHEDULE(static,1) do sh=1,shortcut(0,1) do sh2=sh,shortcut(0,1) exa = 0 @@ -1094,7 +1067,7 @@ subroutine H_S2_u_0_mrcc_nstates(v_0,s_0,u_0,H_jj,S2_jj,n,keys_tmp,Nint,istate_i enddo enddo !$OMP END DO - !$OMP DO SCHEDULE(guided) + !$OMP DO SCHEDULE(static,1) do sh=1,shortcut(0,2) do i=shortcut(sh,2),shortcut(sh+1,2)-1 org_i = sort_idx(i,2) diff --git a/plugins/MRCC_Utils/mrcc_utils.irp.f b/plugins/MRCC_Utils/mrcc_utils.irp.f index e0da2f20..79d139cf 100644 --- a/plugins/MRCC_Utils/mrcc_utils.irp.f +++ b/plugins/MRCC_Utils/mrcc_utils.irp.f @@ -771,10 +771,8 @@ END_PROVIDER factor = 1.d0 resold = huge(1.d0) - do k=0,10*hh_nex + do k=0,hh_nex/4 res = 0.d0 - !$OMP PARALLEL default(shared) private(cx, i, a_col, a_coll) reduction(+:res) - !$OMP DO do a_coll = 1, n_exc_active a_col = active_pp_idx(a_coll) cx = 0.d0 @@ -785,21 +783,20 @@ END_PROVIDER res = res + (X_new(a_col) - X(a_col))*(X_new(a_col) - X(a_col)) X(a_col) = X_new(a_col) end do - !$OMP END DO - !$OMP END PARALLEL if (res > resold) then factor = factor * 0.5d0 endif resold = res - if(iand(k, 4095) == 0) then + if(iand(k, 127) == 0) then print *, "res ", k, res end if if(res < 1d-10) exit end do dIj_unique(1:size(X), s) = X(1:size(X)) + print *, "res ", k, res enddo @@ -831,21 +828,23 @@ END_PROVIDER do s=1,N_states norm = 0.d0 - double precision :: f + double precision :: f, g, gmax + gmax = 1.d0*maxval(dabs(psi_non_ref_coef(:,s))) do i=1,N_det_non_ref - if (rho_mrcc(i,s) == 0.d0) then - rho_mrcc(i,s) = 1.d-32 - endif - if (lambda_type == 2) then f = 1.d0 else + if (rho_mrcc(i,s) == 0.d0) then + cycle + endif ! f is such that f.\tilde{c_i} = c_i f = psi_non_ref_coef(i,s) / rho_mrcc(i,s) ! Avoid numerical instabilities - f = min(f,2.d0) - f = max(f,-2.d0) +! g = 1.d0+dabs(gmax / psi_non_ref_coef(i,s) ) + g = 2.d0+100.d0*exp(-20.d0*dabs(psi_non_ref_coef(i,s)/gmax)) + f = min(f, g) + f = max(f,-g) endif norm = norm + f*f *rho_mrcc(i,s)*rho_mrcc(i,s) @@ -1087,6 +1086,22 @@ end function end do hh_shortcut(hh_shortcut(0)+1) = s+1 + if (hh_shortcut(0) > N_hh_exists) then + print *, 'Error in ', irp_here + print *, 'hh_shortcut(0) :', hh_shortcut(0) + print *, 'N_hh_exists : ', N_hh_exists + print *, 'Is your active space defined?' + stop + endif + + if (hh_shortcut(hh_shortcut(0)+1)-1 > N_pp_exists) then + print *, 'Error 1 in ', irp_here + print *, 'hh_shortcut(hh_shortcut(0)+1)-1 :', hh_shortcut(hh_shortcut(0)+1)-1 + print *, 'N_pp_exists : ', N_pp_exists + print *, 'Is your active space defined?' + stop + endif + do s=2,4,2 do i=1,hh_shortcut(0) if(hh_exists(s, i) == 0) then @@ -1097,6 +1112,7 @@ end function end if end do + do i=1,hh_shortcut(hh_shortcut(0)+1)-1 if(pp_exists(s, i) == 0) then pp_exists(s-1, i) = 0 diff --git a/plugins/mrcepa0/dressing.irp.f b/plugins/mrcepa0/dressing.irp.f index c772e2aa..cebe0a44 100644 --- a/plugins/mrcepa0/dressing.irp.f +++ b/plugins/mrcepa0/dressing.irp.f @@ -23,7 +23,7 @@ use bitmasks !$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_det_non_ref, N_det_ref, delta_ii_mrcc, delta_ij_mrcc, delta_ii_s2_mrcc, delta_ij_s2_mrcc) & - !$OMP private(h, n, mask, omask, buf, ok, iproc) + !$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 @@ -232,12 +232,6 @@ subroutine mrcc_part_dress(delta_ij_, delta_ii_,delta_ij_s2_, delta_ii_s2_,i_gen !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) = 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) @@ -250,6 +244,10 @@ subroutine mrcc_part_dress(delta_ij_, delta_ii_,delta_ij_s2_, delta_ii_s2_,i_gen call apply_excitation(psi_ref(1,1,i_I), exc, tmp_det, ok, Nint) if(.not. ok) cycle + do i_state=1,N_states + dIK(i_state) = dij(i_I, idx_alpha(k_sd), i_state) + enddo + ! do i_state=1,N_states dka(i_state) = 0.d0 @@ -268,12 +266,8 @@ subroutine mrcc_part_dress(delta_ij_, delta_ii_,delta_ij_s2_, delta_ii_s2_,i_gen loop = .false. if (.not.loop) then call get_excitation(psi_ref(1,1,i_I),psi_non_ref(1,1,idx_alpha(l_sd)),exc,degree,phase2,Nint) - 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) = 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 endif @@ -292,7 +286,6 @@ subroutine mrcc_part_dress(delta_ij_, delta_ii_,delta_ij_s2_, delta_ii_s2_,i_gen k_sd = idx_alpha(l_sd) hla = hij_cache(k_sd) sla = sij_cache(k_sd) -! call i_h_j(tq(1,1,i_alpha),psi_non_ref(1,1,idx_alpha(l_sd)),Nint,hla) do i_state=1,N_states dIa_hla(i_state,k_sd) = dIa(i_state) * hla dIa_sla(i_state,k_sd) = dIa(i_state) * sla @@ -336,6 +329,7 @@ end integer :: i, j, i_state !mrmode : 1=mrcepa0, 2=mrsc2 add, 3=mrcc + PROVIDE psi_ref_lock if(mrmode == 3) then do i = 1, N_det_ref diff --git a/plugins/mrcepa0/mrcc.irp.f b/plugins/mrcepa0/mrcc.irp.f index a5614942..bb184761 100644 --- a/plugins/mrcepa0/mrcc.irp.f +++ b/plugins/mrcepa0/mrcc.irp.f @@ -5,7 +5,7 @@ program mrsc2sub !mrmode : 1=mrcepa0, 2=mrsc2 add, 3=mrcc mrmode = 3 - + read_wf = .True. SOFT_TOUCH read_wf call set_generators_bitmasks_as_holes_and_particles