diff --git a/plugins/MRCC_Utils/davidson.irp.f b/plugins/MRCC_Utils/davidson.irp.f index a67ca676..69faf00c 100644 --- a/plugins/MRCC_Utils/davidson.irp.f +++ b/plugins/MRCC_Utils/davidson.irp.f @@ -94,7 +94,6 @@ subroutine davidson_diag_hjj_mrcc(dets_in,u_in,H_jj,energies,dim_in,sze,N_st,N_s double precision, allocatable :: overlap(:,:) double precision :: u_dot_v, u_dot_u - integer, allocatable :: kl_pairs(:,:) integer :: k_pairs, kl integer :: iter2 @@ -144,7 +143,6 @@ subroutine davidson_diag_hjj_mrcc(dets_in,u_in,H_jj,energies,dim_in,sze,N_st,N_s sze_8 = align_double(sze) allocate( & - kl_pairs(2,N_st_diag*(N_st_diag+1)/2), & W(sze_8,N_st_diag,davidson_sze_max), & U(sze_8,N_st_diag,davidson_sze_max), & R(sze_8,N_st_diag), & @@ -360,7 +358,6 @@ subroutine davidson_diag_hjj_mrcc(dets_in,u_in,H_jj,energies,dim_in,sze,N_st,N_s call write_time(iunit) deallocate ( & - kl_pairs, & W, residual_norm, & U, overlap, & R, c, & @@ -649,7 +646,6 @@ subroutine davidson_diag_hjj_sjj_mrcc(dets_in,u_in,H_jj,S2_jj,energies,dim_in,sz double precision, allocatable :: overlap(:,:) double precision :: u_dot_v, u_dot_u - integer, allocatable :: kl_pairs(:,:) integer :: k_pairs, kl integer :: iter2 @@ -661,7 +657,7 @@ subroutine davidson_diag_hjj_sjj_mrcc(dets_in,u_in,H_jj,S2_jj,energies,dim_in,sz character*(16384) :: write_buffer double precision :: to_print(3,N_st) double precision :: cpu, wall - integer :: shift, shift2 + integer :: shift, shift2, itermax include 'constants.include.F' !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: U, W, R, S, y, h, lambda @@ -710,23 +706,30 @@ subroutine davidson_diag_hjj_sjj_mrcc(dets_in,u_in,H_jj,S2_jj,energies,dim_in,sz else delta = 0.d0 endif - - allocate( & - kl_pairs(2,N_st_diag*(N_st_diag+1)/2), & - W(sze_8,N_st_diag*davidson_sze_max), & - U(sze_8,N_st_diag*davidson_sze_max), & - R(sze_8,N_st_diag), & - S(sze_8,N_st_diag*davidson_sze_max), & - h(N_st_diag*davidson_sze_max,N_st_diag*davidson_sze_max), & - y(N_st_diag*davidson_sze_max,N_st_diag*davidson_sze_max), & - s_(N_st_diag*davidson_sze_max,N_st_diag*davidson_sze_max), & - s_tmp(N_st_diag*davidson_sze_max,N_st_diag*davidson_sze_max), & - residual_norm(N_st_diag), & - overlap(N_st_diag,N_st_diag), & - c(N_st_diag*davidson_sze_max), & - s2(N_st_diag*davidson_sze_max), & - lambda(N_st_diag*davidson_sze_max)) + itermax = min(davidson_sze_max, sze/N_st_diag) + allocate( & + W(sze_8,N_st_diag*itermax), & + U(sze_8,N_st_diag*itermax), & + S(sze_8,N_st_diag*itermax), & + h(N_st_diag*itermax,N_st_diag*itermax), & + 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), & + c(N_st_diag*itermax), & + s2(N_st_diag*itermax), & + lambda(N_st_diag*itermax)) + + h = 0.d0 + s_ = 0.d0 + s_tmp = 0.d0 + U = 0.d0 + W = 0.d0 + S = 0.d0 + y = 0.d0 + + ASSERT (N_st > 0) ASSERT (N_st_diag >= N_st) ASSERT (sze > 0) @@ -738,25 +741,25 @@ subroutine davidson_diag_hjj_sjj_mrcc(dets_in,u_in,H_jj,S2_jj,energies,dim_in,sz converged = .False. - do k=1,N_st - call normalize(u_in(1,k),sze) - enddo - - do k=N_st+1,N_st_diag + double precision :: r1, r2 + do k=N_st+1,N_st_diag-2,2 do i=1,sze - double precision :: r1, r2 call random_number(r1) call random_number(r2) - u_in(i,k) = dsqrt(-2.d0*dlog(r1))*dcos(dtwo_pi*r2) + r1 = dsqrt(-2.d0*dlog(r1)) + r2 = dtwo_pi*r2 + u_in(i,k) = r1*dcos(r2) + u_in(i,k+1) = r1*dsin(r2) + enddo + enddo + do k=N_st_diag-1,N_st_diag + do i=1,sze + call random_number(r1) + call random_number(r2) + r1 = dsqrt(-2.d0*dlog(r1)) + r2 = dtwo_pi*r2 + u_in(i,k) = r1*dcos(r2) enddo - - ! Gram-Schmidt - ! ------------ - call dgemv('T',sze,k-1,1.d0,u_in,size(u_in,1), & - u_in(1,k),1,0.d0,c,1) - call dgemv('N',sze,k-1,-1.d0,u_in,size(u_in,1), & - c,1,1.d0,u_in(1,k),1) - call normalize(u_in(1,k),sze) enddo @@ -773,10 +776,10 @@ subroutine davidson_diag_hjj_sjj_mrcc(dets_in,u_in,H_jj,S2_jj,energies,dim_in,sz shift = N_st_diag*(iter-1) shift2 = N_st_diag*iter + call ortho_qr(U,size(U,1),sze,shift2) ! Compute |W_k> = \sum_i |i> ! ----------------------------------------- - call H_S2_u_0_mrcc_nstates(W(1,shift+1),S(1,shift+1),U(1,shift+1),H_jj,S2_jj,sze,dets_in,Nint,& istate,N_st_diag,sze_8) @@ -786,19 +789,6 @@ subroutine davidson_diag_hjj_sjj_mrcc(dets_in,u_in,H_jj,S2_jj,energies,dim_in,sz ! ------------------------------------------- -! do l=1,N_st_diag -! do k=1,N_st_diag -! do iter2=1,iter-1 -! h(k,iter2,l,iter) = u_dot_v(U(1,k,iter2),W(1,l,iter),sze) -! h(k,iter,l,iter2) = h(k,iter2,l,iter) -! enddo -! enddo -! do k=1,l -! h(k,iter,l,iter) = u_dot_v(U(1,k,iter),W(1,l,iter),sze) -! h(l,iter,k,iter) = h(k,iter,l,iter) -! enddo -! enddo - call dgemm('T','N', shift2, N_st_diag, sze, & 1.d0, U, size(U,1), W(1,shift+1), size(W,1), & 0.d0, h(1,shift+1), size(h,1)) @@ -829,7 +819,7 @@ subroutine davidson_diag_hjj_sjj_mrcc(dets_in,u_in,H_jj,S2_jj,energies,dim_in,sz if (s2_eig) then logical :: state_ok(N_st_diag*davidson_sze_max) do k=1,shift2 - state_ok(k) = (dabs(s2(k)-expected_s2) < 0.3d0) + state_ok(k) = (dabs(s2(k)-expected_s2) < 0.6d0) enddo do k=1,shift2 if (.not. state_ok(k)) then @@ -851,22 +841,6 @@ subroutine davidson_diag_hjj_sjj_mrcc(dets_in,u_in,H_jj,S2_jj,energies,dim_in,sz ! Express eigenvectors of h in the determinant basis ! -------------------------------------------------- -! do k=1,N_st_diag -! do i=1,sze -! U(i,shift2+k) = 0.d0 -! W(i,shift2+k) = 0.d0 -! S(i,shift2+k) = 0.d0 -! enddo -! do l=1,N_st_diag*iter -! do i=1,sze -! U(i,shift2+k) = U(i,shift2+k) + U(i,l)*y(l,k) -! W(i,shift2+k) = W(i,shift2+k) + W(i,l)*y(l,k) -! S(i,shift2+k) = S(i,shift2+k) + S(i,l)*y(l,k) -! enddo -! enddo -! enddo -! -! call dgemm('N','N', sze, N_st_diag, shift2, & 1.d0, U, size(U,1), y, size(y,1), 0.d0, U(1,shift2+1), size(U,1)) call dgemm('N','N', sze, N_st_diag, shift2, & @@ -876,83 +850,39 @@ subroutine davidson_diag_hjj_sjj_mrcc(dets_in,u_in,H_jj,S2_jj,energies,dim_in,sz ! Compute residual vector ! ----------------------- - -! do k=1,N_st_diag -! print *, s2(k) -! s2(k) = u_dot_v(U(1,shift2+k), S(1,shift2+k), sze) + S_z2_Sz -! print *, s2(k) -! print *, '' -! pause -! enddo + do k=1,N_st_diag do i=1,sze - R(i,k) = (lambda(k) * U(i,shift2+k) - W(i,shift2+k) ) & - * (1.d0 + s2(k) * U(i,shift2+k) - S(i,shift2+k) - S_z2_Sz) + U(i,shift2+k) = (lambda(k) * U(i,shift2+k) - W(i,shift2+k) ) & + * (1.d0 + s2(k) * U(i,shift2+k) - S(i,shift2+k) - S_z2_Sz & + )/max(H_jj(i) - lambda (k),1.d-2) enddo if (k <= N_st) then - residual_norm(k) = u_dot_u(R(1,k),sze) + residual_norm(k) = u_dot_u(U(1,shift2+k),sze) to_print(1,k) = lambda(k) + nuclear_repulsion to_print(2,k) = s2(k) to_print(3,k) = residual_norm(k) - if (residual_norm(k) > 1.e9) then - stop 'Davidson failed' - endif endif enddo - write(iunit,'(X,I3,X,100(X,F16.10,X,F11.6,X,E11.3))') iter, to_print(:,1:N_st) + write(iunit,'(X,I3,X,100(X,F16.10,X,F11.6,X,E11.3))') iter, to_print(1:3,1:N_st) call davidson_converged(lambda,residual_norm,wall,iter,cpu,N_st,converged) + do k=1,N_st + if (residual_norm(k) > 1.e8) then + print *, '' + stop 'Davidson failed' + endif + enddo if (converged) then exit endif - - ! Davidson step - ! ------------- - - do k=1,N_st_diag - do i=1,sze - U(i,shift2+k) = - R(i,k)/max(H_jj(i) - lambda(k),1.d-2) - enddo - enddo - - ! Gram-Schmidt - ! ------------ - - do k=1,N_st_diag - -! do l=1,N_st_diag*iter -! c(1) = u_dot_v(U(1,shift2+k),U(1,l),sze) -! do i=1,sze -! U(i,k,iter+1) = U(i,shift2+k) - c(1) * U(i,l) -! enddo -! enddo -! - call dgemv('T',sze,N_st_diag*iter,1.d0,U,size(U,1), & - U(1,shift2+k),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,shift2+k),1) -! -! do l=1,k-1 -! c(1) = u_dot_v(U(1,shift2+k),U(1,shift2+l),sze) -! do i=1,sze -! U(i,k,iter+1) = U(i,shift2+k) - c(1) * U(i,shift2+l) -! enddo -! enddo -! - call dgemv('T',sze,k-1,1.d0,U(1,shift2+1),size(U,1), & - U(1,shift2+k),1,0.d0,c,1) - call dgemv('N',sze,k-1,-1.d0,U(1,shift2+1),size(U,1), & - c,1,1.d0,U(1,shift2+k),1) - - call normalize( U(1,shift2+k), sze ) - enddo enddo if (.not.converged) then - iter = davidson_sze_max-1 + iter = itermax-1 endif - + ! Re-contract to u_in ! ----------- @@ -960,15 +890,6 @@ subroutine davidson_diag_hjj_sjj_mrcc(dets_in,u_in,H_jj,S2_jj,energies,dim_in,sz energies(k) = lambda(k) enddo -! do k=1,N_st_diag -! do i=1,sze -! do l=1,iter*N_st_diag -! u_in(i,k) += U(i,l)*y(l,k) -! 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), 0.d0, u_in, size(u_in,1)) @@ -983,7 +904,6 @@ subroutine davidson_diag_hjj_sjj_mrcc(dets_in,u_in,H_jj,S2_jj,energies,dim_in,sz call write_time(iunit) deallocate ( & - kl_pairs, & W, residual_norm, & U, overlap, & R, c, S, & diff --git a/plugins/MRPT_Utils/give_2h2p.irp.f b/plugins/MRPT_Utils/give_2h2p.irp.f new file mode 100644 index 00000000..df71e594 --- /dev/null +++ b/plugins/MRPT_Utils/give_2h2p.irp.f @@ -0,0 +1,35 @@ +subroutine give_2h2p(contrib_2h2p) + implicit none + double precision, intent(out) :: contrib_2h2p(N_states) + integer :: i,j,k,l,m + integer :: iorb,jorb,korb,lorb + + double precision :: get_mo_bielec_integral + double precision :: direct_int,exchange_int + double precision :: numerator,denominator(N_states) + + contrib_2h2p = 0.d0 + do i = 1, n_inact_orb + iorb = list_inact(i) + do j = 1, n_inact_orb + jorb = list_inact(j) + do k = 1, n_virt_orb + korb = list_virt(k) + do l = 1, n_virt_orb + lorb = list_virt(l) + direct_int = get_mo_bielec_integral(iorb,jorb,korb,lorb ,mo_integrals_map) + exchange_int = get_mo_bielec_integral(iorb,jorb,lorb,korb ,mo_integrals_map) + numerator = 3.d0 * direct_int*direct_int + exchange_int*exchange_int -2.d0 * exchange_int * direct_int + do m = 1, N_states + denominator(m) = fock_core_inactive_total_spin_trace(iorb,m) + fock_core_inactive_total_spin_trace(jorb,m) & + -fock_virt_total_spin_trace(korb,m) - fock_virt_total_spin_trace(lorb,m) + contrib_2h2p(m) += numerator / denominator(m) + enddo + enddo + enddo + enddo + enddo + contrib_2h2p = contrib_2h2p*0.5d0 + +end + diff --git a/plugins/MRPT_Utils/mrpt_utils.irp.f b/plugins/MRPT_Utils/mrpt_utils.irp.f index 80739aa2..d7b1f0f6 100644 --- a/plugins/MRPT_Utils/mrpt_utils.irp.f +++ b/plugins/MRPT_Utils/mrpt_utils.irp.f @@ -262,89 +262,87 @@ END_PROVIDER print*, 'Davidson not yet implemented for the dressing ... ' stop - else if (diag_algorithm == "Lapack") then - - allocate (eigenvectors(size(H_matrix_all_dets,1),N_det)) - allocate (eigenvalues(N_det)) - call lapack_diag(eigenvalues,eigenvectors, & - Hmatrix_dressed_pt2_new_symmetrized(1,1,1),N_det,N_det) - CI_electronic_dressed_pt2_new_energy(:) = 0.d0 - if (s2_eig) then - i_state = 0 - allocate (s2_eigvalues(N_det)) - allocate(index_good_state_array(N_det),good_state_array(N_det)) - good_state_array = .False. - do j=1,N_det - call get_s2_u0(psi_det,eigenvectors(1,j),N_det,size(eigenvectors,1),s2) - s2_eigvalues(j) = s2 - ! Select at least n_states states with S^2 values closed to "expected_s2" - if(dabs(s2-expected_s2).le.0.3d0)then - i_state +=1 - index_good_state_array(i_state) = j - good_state_array(j) = .True. - endif - if(i_state.eq.N_states) then - exit - endif - enddo - if(i_state .ne.0)then - ! Fill the first "i_state" states that have a correct S^2 value - do j = 1, i_state - do i=1,N_det - CI_dressed_pt2_new_eigenvectors(i,j) = eigenvectors(i,index_good_state_array(j)) - enddo - CI_electronic_dressed_pt2_new_energy(j) = eigenvalues(index_good_state_array(j)) - CI_dressed_pt2_new_eigenvectors_s2(j) = s2_eigvalues(index_good_state_array(j)) - enddo - i_other_state = 0 - do j = 1, N_det - if(good_state_array(j))cycle - i_other_state +=1 - if(i_state+i_other_state.gt.n_states_diag)then - exit - endif - call get_s2_u0(psi_det,eigenvectors(1,j),N_det,size(eigenvectors,1),s2) - do i=1,N_det - CI_dressed_pt2_new_eigenvectors(i,i_state+i_other_state) = eigenvectors(i,j) - enddo - CI_electronic_dressed_pt2_new_energy(i_state+i_other_state) = eigenvalues(j) - CI_dressed_pt2_new_eigenvectors_s2(i_state+i_other_state) = s2 - enddo - - deallocate(index_good_state_array,good_state_array) + else if (diag_algorithm == "Lapack") then - else - print*,'' - print*,'!!!!!!!! WARNING !!!!!!!!!' - print*,' Within the ',N_det,'determinants selected' - print*,' and the ',N_states_diag,'states requested' - print*,' We did not find any state with S^2 values close to ',expected_s2 - print*,' We will then set the first N_states eigenvectors of the H matrix' - print*,' as the CI_dressed_pt2_new_eigenvectors' - print*,' You should consider more states and maybe ask for diagonalize_s2 to be .True. or just enlarge the CI space' - print*,'' - do j=1,min(N_states_diag,N_det) - do i=1,N_det - CI_dressed_pt2_new_eigenvectors(i,j) = eigenvectors(i,j) - enddo - CI_electronic_dressed_pt2_new_energy(j) = eigenvalues(j) - CI_dressed_pt2_new_eigenvectors_s2(j) = s2_eigvalues(j) + allocate (eigenvectors(size(H_matrix_all_dets,1),N_det)) + allocate (eigenvalues(N_det)) + call lapack_diag(eigenvalues,eigenvectors, & + H_matrix_all_dets,size(H_matrix_all_dets,1),N_det) + CI_electronic_energy(:) = 0.d0 + if (s2_eig) then + i_state = 0 + allocate (s2_eigvalues(N_det)) + allocate(index_good_state_array(N_det),good_state_array(N_det)) + good_state_array = .False. + call u_0_S2_u_0(s2_eigvalues,eigenvectors,N_det,psi_det,N_int,& + N_det,size(eigenvectors,1)) + do j=1,N_det + ! Select at least n_states states with S^2 values closed to "expected_s2" + if(dabs(s2_eigvalues(j)-expected_s2).le.0.5d0)then + i_state +=1 + index_good_state_array(i_state) = j + good_state_array(j) = .True. + endif + if(i_state.eq.N_states) then + exit + endif enddo - endif - deallocate(s2_eigvalues) - else - ! Select the "N_states_diag" states of lowest energy - do j=1,min(N_det,N_states_diag) - call get_s2_u0(psi_det,eigenvectors(1,j),N_det,N_det,s2) - do i=1,N_det - CI_dressed_pt2_new_eigenvectors(i,j) = eigenvectors(i,j) - enddo - CI_electronic_dressed_pt2_new_energy(j) = eigenvalues(j) - CI_dressed_pt2_new_eigenvectors_s2(j) = s2 - enddo - endif - deallocate(eigenvectors,eigenvalues) - endif + if(i_state .ne.0)then + ! Fill the first "i_state" states that have a correct S^2 value + do j = 1, i_state + do i=1,N_det + CI_eigenvectors(i,j) = eigenvectors(i,index_good_state_array(j)) + enddo + CI_electronic_energy(j) = eigenvalues(index_good_state_array(j)) + CI_eigenvectors_s2(j) = s2_eigvalues(index_good_state_array(j)) + enddo + i_other_state = 0 + do j = 1, N_det + if(good_state_array(j))cycle + i_other_state +=1 + if(i_state+i_other_state.gt.n_states_diag)then + exit + endif + do i=1,N_det + CI_eigenvectors(i,i_state+i_other_state) = eigenvectors(i,j) + enddo + CI_electronic_energy(i_state+i_other_state) = eigenvalues(j) + CI_eigenvectors_s2(i_state+i_other_state) = s2_eigvalues(i_state+i_other_state) + enddo + + else + print*,'' + print*,'!!!!!!!! WARNING !!!!!!!!!' + print*,' Within the ',N_det,'determinants selected' + print*,' and the ',N_states_diag,'states requested' + print*,' We did not find any state with S^2 values close to ',expected_s2 + print*,' We will then set the first N_states eigenvectors of the H matrix' + print*,' as the CI_eigenvectors' + print*,' You should consider more states and maybe ask for s2_eig to be .True. or just enlarge the CI space' + print*,'' + do j=1,min(N_states_diag,N_det) + do i=1,N_det + CI_eigenvectors(i,j) = eigenvectors(i,j) + enddo + CI_electronic_energy(j) = eigenvalues(j) + CI_eigenvectors_s2(j) = s2_eigvalues(j) + enddo + endif + deallocate(index_good_state_array,good_state_array) + deallocate(s2_eigvalues) + else + call u_0_S2_u_0(CI_eigenvectors_s2,eigenvectors,N_det,psi_det,N_int,& + min(N_det,N_states_diag),size(eigenvectors,1)) + ! Select the "N_states_diag" states of lowest energy + do j=1,min(N_det,N_states_diag) + do i=1,N_det + CI_eigenvectors(i,j) = eigenvectors(i,j) + enddo + CI_electronic_energy(j) = eigenvalues(j) + enddo + endif + deallocate(eigenvectors,eigenvalues) + endif END_PROVIDER