diff --git a/plugins/DavidsonDressed/diagonalize_CI.irp.f b/plugins/DavidsonDressed/diagonalize_CI.irp.f index 2ee540e2..ddea2950 100644 --- a/plugins/DavidsonDressed/diagonalize_CI.irp.f +++ b/plugins/DavidsonDressed/diagonalize_CI.irp.f @@ -63,8 +63,8 @@ END_PROVIDER call davidson_diag_HS2(psi_det,CI_eigenvectors_dressed, CI_eigenvectors_s2_dressed,& size(CI_eigenvectors_dressed,1), CI_electronic_energy_dressed,& N_det,min(N_det,N_states),min(N_det,N_states_diag),N_int,1) - call u_0_S2_u_0(CI_eigenvectors_s2_dressed,CI_eigenvectors_dressed,N_det,psi_det,N_int,& - N_states_diag,size(CI_eigenvectors_dressed,1)) +! call u_0_S2_u_0(CI_eigenvectors_s2_dressed,CI_eigenvectors_dressed,N_det,psi_det,N_int,& +! N_states_diag,size(CI_eigenvectors_dressed,1)) else if (diag_algorithm == "Lapack") then diff --git a/plugins/dress_zmq/dress_general.irp.f b/plugins/dress_zmq/dress_general.irp.f index 068d811e..731011e3 100644 --- a/plugins/dress_zmq/dress_general.irp.f +++ b/plugins/dress_zmq/dress_general.irp.f @@ -29,27 +29,26 @@ subroutine run_dressing(N_st,energy) E_new = 0.d0 delta_E = 1.d0 iteration = 0 - do while (delta_E > thresh_dress) - iteration += 1 + do iteration=1,n_it_dress_max print *, '===============================================' print *, 'Iteration', iteration, '/', n_it_dress_max print *, '===============================================' print *, '' - E_old = dress_e0_denominator(1) !sum(ci_energy_dressed(1:N_states)) + E_old = sum(psi_energy(:)) do i=1,N_st call write_double(6,ci_energy_dressed(i),"Energy") enddo call diagonalize_ci_dressed - E_new = dress_e0_denominator(1) !sum(ci_energy_dressed(1:N_states)) + E_new = sum(psi_energy(:)) delta_E = (E_new - E_old)/dble(N_states) print *, '' call write_double(6,thresh_dress,"thresh_dress") - call write_double(6,delta_E,"delta_E") + call write_double(6,delta_E,"delta_E (undressed)") delta_E = dabs(delta_E) call save_wavefunction ! call ezfio_set_dress_zmq_energy(ci_energy_dressed(1)) - if (iteration >= n_it_dress_max) then + if (delta_E < thresh_dress) then exit endif enddo diff --git a/plugins/dress_zmq/dressing.irp.f b/plugins/dress_zmq/dressing.irp.f index 0c15ee0b..928c72c5 100644 --- a/plugins/dress_zmq/dressing.irp.f +++ b/plugins/dress_zmq/dressing.irp.f @@ -78,7 +78,7 @@ BEGIN_PROVIDER [ double precision, delta_ij, (N_states,N_det,2) ] delta_ij = 0d0 - E_CI_before(:) = dress_E0_denominator(:) + nuclear_repulsion + E_CI_before(:) = psi_energy(:) + nuclear_repulsion threshold_selectors = 1.d0 threshold_generators = 1d0 ! if(errr /= 0d0) then diff --git a/plugins/dress_zmq/dressing_vector.irp.f b/plugins/dress_zmq/dressing_vector.irp.f index 5a8fee3b..e098cf78 100644 --- a/plugins/dress_zmq/dressing_vector.irp.f +++ b/plugins/dress_zmq/dressing_vector.irp.f @@ -18,12 +18,6 @@ dressing_column_h(j,k) = delta_ij(k,j,1) dressing_column_s(j,k) = delta_ij(k,j,2) enddo -! tmp = u_dot_v(dressing_column_h(1,k), psi_coef(1,k), N_det) & -! - dressing_column_h(l,k) * psi_coef(l,k) -! dressing_column_h(l,k) -= tmp * f -! tmp = u_dot_v(dressing_column_s(1,k), psi_coef(1,k), N_det) & -! - dressing_column_s(l,k) * psi_coef(l,k) -! dressing_column_s(l,k) -= tmp * f enddo END_PROVIDER diff --git a/plugins/dress_zmq/energy.irp.f b/plugins/dress_zmq/energy.irp.f index b8948219..7e223a2e 100644 --- a/plugins/dress_zmq/energy.irp.f +++ b/plugins/dress_zmq/energy.irp.f @@ -13,13 +13,24 @@ BEGIN_PROVIDER [ double precision, dress_E0_denominator, (N_states) ] END_DOC integer :: i if (initialize_dress_E0_denominator) then - call u_0_H_u_0(dress_E0_denominator,psi_coef,N_det,psi_det,N_int,N_states,size(psi_coef,1)) - do i=N_det+1,N_states - dress_E0_denominator(i) = 0.d0 - enddo + if (h0_type == "EN") then + dress_E0_denominator(1:N_states) = psi_energy(1:N_states) + else if (h0_type == "Barycentric") then +! dress_E0_denominator(1:N_states) = barycentric_electronic_energy(1:N_states) + dress_E0_denominator(1:N_states) = minval(diagonal_H_matrix_on_psi_det(1:N_det)) + else + print *, h0_type, ' not implemented' + stop + endif +! call u_0_H_u_0(dress_E0_denominator,psi_coef,N_det,psi_det,N_int,N_states,size(psi_coef,1)) +! do i=N_det+1,N_states +! dress_E0_denominator(i) = 0.d0 +! enddo call write_double(6,dress_E0_denominator(1)+nuclear_repulsion, 'dress Energy denominator') else dress_E0_denominator = -huge(1.d0) endif END_PROVIDER + + diff --git a/plugins/shiftedbk/shifted_bk_routines.irp.f b/plugins/shiftedbk/shifted_bk_routines.irp.f index e88d153c..80383e39 100644 --- a/plugins/shiftedbk/shifted_bk_routines.irp.f +++ b/plugins/shiftedbk/shifted_bk_routines.irp.f @@ -4,6 +4,7 @@ &BEGIN_PROVIDER [ double precision, a_s2_i, (N_det, Nproc) ] implicit none current_generator_(:) = 0 + fock_diag_tmp_(:,:,:) = 0.d0 a_h_i = 0d0 a_s2_i = 0d0 END_PROVIDER @@ -52,7 +53,7 @@ subroutine dress_with_alpha_buffer(Nstates,Ndet,Nint,delta_ij_loc, i_gen, minili do i=1,Nstates - de = E0_denominator(i) - haa + de = dress_E0_denominator(i) - haa if(DABS(de) < 1D-5) cycle c_alpha = a_h_psi(i) / de @@ -76,23 +77,4 @@ BEGIN_PROVIDER [ logical, initialize_E0_denominator ] END_PROVIDER -BEGIN_PROVIDER [ double precision, E0_denominator, (N_states) ] - implicit none - BEGIN_DOC - ! E0 in the denominator of the PT2 - END_DOC - if (initialize_E0_denominator) then - if (h0_type == "EN") then - E0_denominator(1:N_states) = psi_energy(1:N_states) - else if (h0_type == "Barycentric") then - E0_denominator(1:N_states) = barycentric_electronic_energy(1:N_states) - else - print *, h0_type, ' not implemented' - stop - endif - else - E0_denominator = -huge(1.d0) - endif -END_PROVIDER - diff --git a/src/Davidson/diagonalization_hs2_dressed.irp.f b/src/Davidson/diagonalization_hs2_dressed.irp.f index 59e9c9fe..c57312a9 100644 --- a/src/Davidson/diagonalization_hs2_dressed.irp.f +++ b/src/Davidson/diagonalization_hs2_dressed.irp.f @@ -35,7 +35,7 @@ subroutine davidson_diag_hs2(dets_in,u_in,s2_out,dim_in,energies,sze,N_st,N_st_d double precision, intent(inout) :: u_in(dim_in,N_st_diag) double precision, intent(out) :: energies(N_st_diag), s2_out(N_st_diag) integer, intent(in) :: dressing_state - double precision, allocatable :: H_jj(:), S2_jj(:) + double precision, allocatable :: H_jj(:) double precision, external :: diag_H_mat_elem, diag_S_mat_elem integer :: i,k @@ -44,7 +44,7 @@ subroutine davidson_diag_hs2(dets_in,u_in,s2_out,dim_in,energies,sze,N_st,N_st_d ASSERT (Nint > 0) ASSERT (Nint == N_int) PROVIDE mo_bielec_integrals_in_map - allocate(H_jj(sze),S2_jj(sze)) + allocate(H_jj(sze)) H_jj(1) = diag_h_mat_elem(dets_in(1,1,1),Nint) !$OMP PARALLEL DEFAULT(NONE) & @@ -60,13 +60,13 @@ subroutine davidson_diag_hs2(dets_in,u_in,s2_out,dim_in,energies,sze,N_st,N_st_d if (dressing_state > 0) then do k=1,N_st do i=1,sze - H_jj(i) += u_in(i,k) * dressing_column_h(i,k) + H_jj(i) += u_in(i,k) * dressing_column_h(i,k) enddo enddo endif call davidson_diag_hjj_sjj(dets_in,u_in,H_jj,S2_out,energies,dim_in,sze,N_st,N_st_diag,Nint,dressing_state) - deallocate (H_jj,S2_jj) + deallocate (H_jj) end @@ -254,9 +254,9 @@ subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,s2_out,energies,dim_in,sze,N_ dressing_column_h, size(dressing_column_h,1), s_tmp, size(s_tmp,1), & 1.d0, W(1,shift+1), size(W,1)) - call dgemm('N','N', sze, N_st_diag, N_st, 0.5d0, & - dressing_column_s, size(dressing_column_s,1), s_tmp, size(s_tmp,1), & - 1.d0, S(1,shift+1), size(S,1)) +! call dgemm('N','N', sze, N_st_diag, N_st, 0.5d0, & +! dressing_column_s, size(dressing_column_s,1), s_tmp, size(s_tmp,1), & +! 1.d0, S(1,shift+1), size(S,1)) call dgemm('T','N', N_st, N_st_diag, sze, 1.d0, & @@ -267,13 +267,13 @@ subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,s2_out,energies,dim_in,sze,N_ psi_coef, size(psi_coef,1), s_tmp, size(s_tmp,1), & 1.d0, W(1,shift+1), size(W,1)) - call dgemm('T','N', N_st, N_st_diag, sze, 1.d0, & - dressing_column_s, size(dressing_column_s,1), & - U(1,shift+1), size(U,1), 0.d0, s_tmp, size(s_tmp,1)) - - call dgemm('N','N', sze, N_st_diag, N_st, 0.5d0, & - psi_coef, size(psi_coef,1), s_tmp, size(s_tmp,1), & - 1.d0, S(1,shift+1), size(S,1)) +! call dgemm('T','N', N_st, N_st_diag, sze, 1.d0, & +! dressing_column_s, size(dressing_column_s,1), & +! U(1,shift+1), size(U,1), 0.d0, s_tmp, size(s_tmp,1)) +! +! call dgemm('N','N', sze, N_st_diag, N_st, 0.5d0, & +! psi_coef, size(psi_coef,1), s_tmp, size(s_tmp,1), & +! 1.d0, S(1,shift+1), size(S,1)) endif diff --git a/src/Determinants/slater_rules.irp.f b/src/Determinants/slater_rules.irp.f index ee597720..a3120be9 100644 --- a/src/Determinants/slater_rules.irp.f +++ b/src/Determinants/slater_rules.irp.f @@ -509,11 +509,12 @@ subroutine i_H_j_s2(key_i,key_j,Nint,hij,s2) select case (degree) case (2) call get_double_excitation(key_i,key_j,exc,phase,Nint) - + ! Mono alpha, mono beta 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 + endif if(exc(1,1,1) == exc(1,2,2) )then - if(exc(1,1,2) == exc(1,2,1)) s2 = -phase !!!!! hij = phase * big_array_exchange_integrals(exc(1,1,1),exc(1,1,2),exc(1,2,1)) else if (exc(1,2,1) ==exc(1,1,2))then hij = phase * big_array_exchange_integrals(exc(1,2,1),exc(1,1,1),exc(1,2,2)) @@ -524,8 +525,8 @@ subroutine i_H_j_s2(key_i,key_j,Nint,hij,s2) exc(1,2,1), & exc(1,2,2) ,mo_integrals_map) endif + ! Double alpha else if (exc(0,1,1) == 2) then - ! Double alpha hij = phase*(get_mo_bielec_integral( & exc(1,1,1), & exc(2,1,1), & @@ -536,8 +537,8 @@ subroutine i_H_j_s2(key_i,key_j,Nint,hij,s2) exc(2,1,1), & exc(2,2,1), & exc(1,2,1) ,mo_integrals_map) ) + ! Double beta else if (exc(0,1,2) == 2) then - ! Double beta hij = phase*(get_mo_bielec_integral( & exc(1,1,2), & exc(2,1,2), & @@ -553,13 +554,13 @@ subroutine i_H_j_s2(key_i,key_j,Nint,hij,s2) call get_mono_excitation(key_i,key_j,exc,phase,Nint) !DIR$ FORCEINLINE call bitstring_to_list_ab(key_i, occ, n_occ_ab, Nint) + ! Mono alpha if (exc(0,1,1) == 1) then - ! Mono alpha m = exc(1,1,1) p = exc(1,2,1) spin = 1 + ! Mono beta else - ! Mono beta m = exc(1,1,2) p = exc(1,2,2) spin = 2 @@ -567,7 +568,7 @@ subroutine i_H_j_s2(key_i,key_j,Nint,hij,s2) call get_mono_excitation_from_fock(key_i,key_j,p,m,spin,phase,hij) case (0) - print *," ZERO" + print *,irp_here,": ZERO" double precision, external :: diag_S_mat_elem s2 = diag_S_mat_elem(key_i,Nint) hij = diag_H_mat_elem(key_i,Nint)