diff --git a/src/ccsd/ccsd_t_space_orb_stoch.irp.f b/src/ccsd/ccsd_t_space_orb_stoch.irp.f index 13fa4f1a..293baa2d 100644 --- a/src/ccsd/ccsd_t_space_orb_stoch.irp.f +++ b/src/ccsd/ccsd_t_space_orb_stoch.irp.f @@ -181,8 +181,8 @@ subroutine ccsd_par_t_space_stoch(nO,nV,t1,t2,f_o,f_v,v_vvvo,v_vvoo,v_vooo,energ integer :: nbuckets nbuckets = 100 + double precision, allocatable :: ED(:) double precision, allocatable :: wsum(:) - allocate(wsum(nbuckets)) converged = .False. Ncomputed = 0_8 @@ -197,7 +197,8 @@ subroutine ccsd_par_t_space_stoch(nO,nV,t1,t2,f_o,f_v,v_vvvo,v_vvoo,v_vooo,energ iright = Nabc integer*8, allocatable :: bounds(:,:) - allocate (bounds(2,nbuckets)) + allocate(wsum(nbuckets), ED(nbuckets), bounds(2,nbuckets)) + ED(:) = 0.d0 do isample=1,nbuckets eta = 1.d0/dble(nbuckets) * dble(isample) ieta = binary_search(waccu,eta,Nabc) @@ -233,7 +234,7 @@ subroutine ccsd_par_t_space_stoch(nO,nV,t1,t2,f_o,f_v,v_vvvo,v_vvoo,v_vooo,energ enddo ! Deterministic part - if (imin < Nabc) then + if (imin <= Nabc) then ieta=imin sampled(ieta) = 0_8 a = abc(1,ieta) @@ -254,7 +255,7 @@ subroutine ccsd_par_t_space_stoch(nO,nV,t1,t2,f_o,f_v,v_vvvo,v_vvoo,v_vooo,energ ! Stochastic part call random_number(eta) do isample=1,nbuckets - if (imin >= bounds(2,isample)) then + if (imin > bounds(2,isample)) then cycle endif ieta = binary_search(waccu,(eta + dble(isample-1))/dble(nbuckets),Nabc)+1 @@ -280,7 +281,7 @@ subroutine ccsd_par_t_space_stoch(nO,nV,t1,t2,f_o,f_v,v_vvvo,v_vvoo,v_vooo,energ enddo call wall_time(t01) - if ((t01-t00 > 1.0d0).or.(imin >= Nabc)) then + if ((t01-t00 > 1.0d0).or.(imin > Nabc)) then !$OMP TASKWAIT call wall_time(t01) @@ -300,8 +301,11 @@ subroutine ccsd_par_t_space_stoch(nO,nV,t1,t2,f_o,f_v,v_vvvo,v_vvoo,v_vooo,energ do isample=1,nbuckets - if (imin >= bounds(2,isample)) then - energy_det = energy_det + sum(memo(bounds(1,isample):bounds(2,isample))) + if (imin > bounds(2,isample)) then + if (ED(isample) == 0.d0) then + ED(isample) = sum(memo(bounds(1,isample):bounds(2,isample))) + endif + energy_det = energy_det + ED(isample) scale = scale - wsum(isample) else exit @@ -310,12 +314,14 @@ subroutine ccsd_par_t_space_stoch(nO,nV,t1,t2,f_o,f_v,v_vvvo,v_vvoo,v_vooo,energ isample = min(isample,nbuckets) do ieta=bounds(1,isample), Nabc - w = dble(max(sampled(ieta),0_8)) - tmp = w * memo(ieta) * Pabc(ieta) - ET = ET + tmp - ET2 = ET2 + tmp * memo(ieta) * Pabc(ieta) - norm = norm + w + if (sampled(ieta) < 0_8) cycle + w = dble(sampled(ieta)) + tmp = w * memo(ieta) * Pabc(ieta) + ET = ET + tmp + ET2 = ET2 + tmp * memo(ieta) * Pabc(ieta) + norm = norm + w enddo + norm = norm/scale if (norm > 0.d0) then energy_stoch = ET / norm @@ -327,7 +333,7 @@ subroutine ccsd_par_t_space_stoch(nO,nV,t1,t2,f_o,f_v,v_vvvo,v_vvoo,v_vooo,energ print '('' '',F20.8, '' '', ES12.4,'' '', F8.2,'' '')', eccsd+energy, dsqrt(variance/(norm-1.d0)), 100.*real(Ncomputed)/real(Nabc) endif !$OMP END MASTER - if (imin >= Nabc) exit + if (imin > Nabc) exit enddo !$OMP END PARALLEL diff --git a/src/davidson/diagonalization_hs2_dressed.irp.f b/src/davidson/diagonalization_hs2_dressed.irp.f index 3513f215..fd967ecc 100644 --- a/src/davidson/diagonalization_hs2_dressed.irp.f +++ b/src/davidson/diagonalization_hs2_dressed.irp.f @@ -139,7 +139,7 @@ subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,s2_out,energies,dim_in,sze,N_ integer :: iter2, itertot double precision, allocatable :: y(:,:), h(:,:), h_p(:,:), lambda(:), s2(:) real, allocatable :: y_s(:,:) - double precision, allocatable :: s_(:,:), s_tmp(:,:) + double precision, allocatable :: s_(:,:), s_tmp(:,:), prev_y(:,:) double precision :: diag_h_mat_elem double precision, allocatable :: residual_norm(:) character*(16384) :: write_buffer @@ -288,6 +288,7 @@ subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,s2_out,energies,dim_in,sze,N_ h(N_st_diag*itermax,N_st_diag*itermax), & ! h_p(N_st_diag*itermax,N_st_diag*itermax), & y(N_st_diag*itermax,N_st_diag*itermax), & + prev_y(N_st_diag*itermax,N_st_diag*itermax), & s_(N_st_diag*itermax,N_st_diag*itermax), & s_tmp(N_st_diag*itermax,N_st_diag*itermax), & residual_norm(N_st_diag), & @@ -301,6 +302,10 @@ subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,s2_out,energies,dim_in,sze,N_ s_ = 0.d0 s_tmp = 0.d0 + prev_y = 0.d0 + do i = 1, N_st_diag*itermax + prev_y(i,i) = 1d0 + enddo ASSERT (N_st > 0) ASSERT (N_st_diag >= N_st) @@ -479,6 +484,10 @@ subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,s2_out,energies,dim_in,sze,N_ if (info > 0) then ! Numerical errors propagate. We need to reduce the number of iterations itermax = iter-1 + + ! eigenvectors of the previous iteration + y = prev_y + shift2 = shift2 - N_st_diag exit endif @@ -553,7 +562,7 @@ subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,s2_out,energies,dim_in,sze,N_ do l = 1, N_st do k = 1, N_st_diag do i = 1, sze - overlp(k+j-1,l) += U(i,l) * U(i,shift2+k) + overlp(k+j-1,l) += u_in(i,l) * U(i,shift2+k) enddo enddo enddo @@ -576,7 +585,7 @@ subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,s2_out,energies,dim_in,sze,N_ enddo ! Maximum overlap - if (dabs(overlp(k,l)) > omax .and. .not. used .and. state_ok(k)) then + if ((dabs(overlp(k,l)) > omax) .and. (.not. used) .and. state_ok(k)) then omax = dabs(overlp(k,l)) idx = k endif @@ -615,6 +624,9 @@ subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,s2_out,energies,dim_in,sze,N_ endif enddo + ! Swapped eigenvectors + prev_y = y + ! if (state_following) then ! ! overlap = -1.d0 @@ -677,7 +689,7 @@ subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,s2_out,energies,dim_in,sze,N_ do i=1,sze U(i,shift2+k) = & (lambda(k) * U(i,shift2+k) - W(i,shift2+k) ) & - /max(H_jj(i) - lambda (k),1.d-2) + /max(dabs(H_jj(i) - lambda (k)),1.d-2) * dsign(1d0,H_jj(i) - lambda (k)) enddo if (k <= N_st) then @@ -792,7 +804,7 @@ subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,s2_out,energies,dim_in,sze,N_ residual_norm, & U, overlap, & h, y_s, S_d, & - y, s_, s_tmp, & + y, s_, s_tmp, prev_y, & lambda & ) FREE nthreads_davidson