From 77f38a94a2924c5f39cdaddd1762e9432e3212bb Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Fri, 14 Apr 2017 11:09:55 +0200 Subject: [PATCH] working on davidson --- src/Davidson/u0Hu0.irp.f | 279 +++++++++++++++--------- src/Determinants/slater_rules.irp.f | 5 +- src/Determinants/spindeterminants.irp.f | 12 + 3 files changed, 195 insertions(+), 101 deletions(-) diff --git a/src/Davidson/u0Hu0.irp.f b/src/Davidson/u0Hu0.irp.f index af01eba8..b7d47e08 100644 --- a/src/Davidson/u0Hu0.irp.f +++ b/src/Davidson/u0Hu0.irp.f @@ -632,9 +632,9 @@ subroutine H_S2_u_0_nstates_new(v_0,s_0,N_st,sze_8) PROVIDE ref_bitmask_energy double precision :: hij, s2 - integer :: i,j + integer :: i,j,k integer :: k_a, k_b, l_a, l_b, m_a, m_b - integer :: degree, istate + integer :: istate integer :: krow, kcol, krow_b, kcol_b integer :: lrow, lcol integer :: mrow, mcol @@ -646,16 +646,20 @@ subroutine H_S2_u_0_nstates_new(v_0,s_0,N_st,sze_8) double precision :: ck(N_st), cl(N_st), cm(N_st) integer :: n_singles, n_doubles integer, allocatable :: singles(:), doubles(:) + integer, allocatable :: singles_a(:,:), singles_b(:,:) integer, allocatable :: idx(:), idx0(:) logical, allocatable :: is_single_a(:) + logical, allocatable :: is_single_b(:) + integer :: maxab, n_singles_max - allocate( buffer(N_int,N_det_alpha_unique), & - singles(N_det_alpha_unique), doubles(N_det_alpha_unique), & + maxab = max(N_det_alpha_unique, N_det_beta_unique) + allocate( buffer(N_int,maxab), & + singles(maxab), doubles(maxab), & is_single_a(N_det_alpha_unique), & - idx(N_det_alpha_unique), idx0(N_det_alpha_unique) ) + is_single_b(N_det_beta_unique), & + idx(maxab), idx0(maxab)) v_0 = 0.d0 - do k_a=1,N_det ! Initial determinant is at k_a in alpha-major representation @@ -690,12 +694,13 @@ subroutine H_S2_u_0_nstates_new(v_0,s_0,N_st,sze_8) i=1 l_a = k_a+1 lcol = psi_bilinear_matrix_columns(l_a) - do while ( (lcol == kcol).and.(l_a <= N_det) ) + do while (lcol == kcol) lrow = psi_bilinear_matrix_rows(l_a) buffer(1:N_int,i) = psi_det_alpha_unique(1:N_int, lrow) idx(i) = lrow - i=i+1 - l_a = l_a + 1 + i = i +1 + l_a = l_a+1 + if (l_a > N_det) exit lcol = psi_bilinear_matrix_columns(l_a) enddo i = i-1 @@ -747,12 +752,13 @@ subroutine H_S2_u_0_nstates_new(v_0,s_0,N_st,sze_8) i=1 l_b = k_b+1 lrow = psi_bilinear_matrix_transp_rows(l_b) - do while ( (lrow == krow).and.(l_b <= N_det) ) + do while (lrow == krow) lcol = psi_bilinear_matrix_transp_columns(l_b) buffer(1:N_int,i) = psi_det_beta_unique(1:N_int, lcol) idx(i) = lcol - i=i+1 - l_b = l_b + 1 + i = i +1 + l_b = l_b+1 + if (l_b > N_det) exit lrow = psi_bilinear_matrix_transp_rows(l_b) enddo i = i-1 @@ -801,115 +807,190 @@ subroutine H_S2_u_0_nstates_new(v_0,s_0,N_st,sze_8) ! Alpha/Beta double excitations ! ============================= - do i=1,N_det_beta_unique + do i=1,maxab idx0(i) = i enddo - is_single_a(:) = .False. - k_a=1 - do i=1,N_det_beta_unique - - ! Select a beta determinant - ! ------------------------- - - spindet(1:N_int) = psi_det_beta_unique(1:N_int, i) - tmp_det(1:N_int,2) = spindet(1:N_int) + ! Prepare the array of all alpha single excitations + ! ------------------------------------------------- + n_singles_max = 0 + do i=1,N_det_alpha_unique + spindet(1:N_int) = psi_det_alpha_unique(1:N_int, i) call get_all_spin_singles( & - psi_det_beta_unique, idx0, spindet, N_int, N_det_beta_unique, & - singles, n_singles ) + psi_det_alpha_unique, idx0, spindet, N_int, N_det_alpha_unique,& + singles, n_singles) + n_singles_max = max(n_singles_max, n_singles) + enddo - do j=1,n_singles - is_single_a( singles(j) ) = .True. - enddo + allocate (singles_a(0:n_singles_max, N_det_alpha_unique)) + do i=1,N_det_alpha_unique + spindet(1:N_int) = psi_det_alpha_unique(1:N_int, i) + call get_all_spin_singles( & + psi_det_alpha_unique, idx0, spindet, N_int, N_det_alpha_unique,& + singles_a(1,i), singles_a(0,i)) + enddo - ! For all alpha.beta pairs with the selected beta - ! ----------------------------------------------- + do k_a=1,N_det + ! Initial determinant is at k_a in alpha-major representation + ! ----------------------------------------------------------------------- + + krow = psi_bilinear_matrix_rows(k_a) kcol = psi_bilinear_matrix_columns(k_a) - do while (kcol < i) - k_a = k_a+1 - if (k_a > N_det) exit - kcol = psi_bilinear_matrix_columns(k_a) + + tmp_det(1:N_int,1) = psi_det_alpha_unique(1:N_int, krow) + tmp_det(1:N_int,2) = psi_det_beta_unique (1:N_int, kcol) + + is_single_a = .False. + do k=1,singles_a(0,krow) + is_single_a( singles_a(k,krow) ) = .True. enddo - do while (kcol == i) + if (k_a > 1) then + if (kcol /= psi_bilinear_matrix_columns(k_a-1)) then + call get_all_spin_singles( & + psi_det_beta_unique, idx0, tmp_det(1,2), N_int, N_det_beta_unique,& + singles, n_singles) + endif + else + call get_all_spin_singles( & + psi_det_beta_unique, idx0, tmp_det(1,2), N_int, N_det_beta_unique,& + singles, n_singles) + endif - krow = psi_bilinear_matrix_rows(k_a) - tmp_det(1:N_int,1) = psi_det_alpha_unique(1:N_int,krow) + ! Loop over singly excited beta columns + ! ------------------------------------- - ! Loop over all alpha.beta pairs with a single exc alpha - ! ------------------------------------------------------ + do i=1,n_singles + lcol = singles(i) + ! TODO cycle if lcol <= kcol + tmp_det2(1:N_int,2) = psi_det_beta_unique(1:N_int, lcol) - l_a = k_a+1 - if (l_a > N_det) exit - lrow = psi_bilinear_matrix_rows(l_a) - lcol = psi_bilinear_matrix_columns(l_a) + l_a = psi_bilinear_matrix_columns_loc(lcol) + ! TODO loop - do while (lrow == krow) - - ! Loop over all alpha.beta pairs with a single exc alpha - ! ------------------------------------------------------ + do while ( l_a < psi_bilinear_matrix_columns_loc(lcol+1) ) + lrow = psi_bilinear_matrix_rows(l_a) if (is_single_a(lrow)) then + tmp_det2(1:N_int,1) = psi_det_alpha_unique(1:N_int, lrow) - tmp_det2(1:N_int,1) = psi_det_alpha_unique(1:N_int,lrow) - - ! Build list of singly excited beta - ! --------------------------------- + call i_H_j_double_alpha_beta(tmp_det,tmp_det2,N_int,hij) + v_0(k_a, 1:N_st) += hij * psi_bilinear_matrix_values(l_a,1:N_st) + endif + l_a += 1 - m_b = psi_bilinear_matrix_order_reverse(l_a) - m_b = m_b+1 - j=1 - do while ( (mrow == lrow) ) - mcol = psi_bilinear_matrix_transp_columns(m_b) - buffer(1:N_int,j) = psi_det_beta_unique(1:N_int,mcol) - idx(j) = mcol - j = j+1 - m_b = m_b+1 - if (m_b <= N_det) exit - mrow = psi_bilinear_matrix_transp_rows(m_b) - enddo - j=j-1 + enddo + enddo - call get_all_spin_singles( & - buffer, idx, tmp_det(1,2), N_int, j, & - doubles, n_doubles) + enddo - ! Compute Hij for all doubles - ! --------------------------- + !---- - m_b = psi_bilinear_matrix_order(l_a)+1 - mcol = psi_bilinear_matrix_transp_columns(m_b) - do j=1,n_doubles - tmp_det2(1:N_int,2) = psi_det_beta_unique(1:N_int, doubles(j) ) - call i_H_j_double_alpha_beta(tmp_det,tmp_det2,N_int,hij) - do while (mcol /= doubles(j)) - m_b = m_b+1 - if (m_b > N_det) exit - mcol = psi_bilinear_matrix_transp_columns(m_b) - enddo - m_a = psi_bilinear_matrix_order_reverse(m_b) +! k_a=1 +! do i=1,N_det_beta_unique +! +! ! Select a beta determinant +! ! ------------------------- +! +! spindet(1:N_int) = psi_det_beta_unique(1:N_int, i) +! tmp_det(1:N_int,2) = spindet(1:N_int) +! +! call get_all_spin_singles( & +! psi_det_beta_unique, idx0, spindet, N_int, N_det_beta_unique, & +! singles, n_singles ) +! +! do j=1,n_singles +! is_single_a( singles(j) ) = .True. +! enddo +! +! ! For all alpha.beta pairs with the selected beta +! ! ----------------------------------------------- +! +! kcol = psi_bilinear_matrix_columns(k_a) +! do while (kcol < i) +! k_a = k_a+1 +! if (k_a > N_det) exit +! kcol = psi_bilinear_matrix_columns(k_a) +! enddo +! +! do while (kcol == i) +! +! krow = psi_bilinear_matrix_rows(k_a) +! tmp_det(1:N_int,1) = psi_det_alpha_unique(1:N_int,krow) +! +! ! Loop over all alpha.beta pairs with a single exc alpha +! ! ------------------------------------------------------ +! +! l_a = k_a+1 +! if (l_a > N_det) exit +! lrow = psi_bilinear_matrix_rows(l_a) +! lcol = psi_bilinear_matrix_columns(l_a) +! +! do while (lrow == krow) +! +! ! Loop over all alpha.beta pairs with a single exc alpha +! ! ------------------------------------------------------ +! if (is_single_a(lrow)) then +! +! tmp_det2(1:N_int,1) = psi_det_alpha_unique(1:N_int,lrow) +! +! ! Build list of singly excited beta +! ! --------------------------------- +! +! m_b = psi_bilinear_matrix_order_reverse(l_a) +! m_b = m_b+1 +! j=1 +! do while ( (mrow == lrow) ) +! mcol = psi_bilinear_matrix_transp_columns(m_b) +! buffer(1:N_int,j) = psi_det_beta_unique(1:N_int,mcol) +! idx(j) = mcol +! j = j+1 +! m_b = m_b+1 +! if (m_b <= N_det) exit +! mrow = psi_bilinear_matrix_transp_rows(m_b) +! enddo +! j=j-1 +! +! call get_all_spin_singles( & +! buffer, idx, tmp_det(1,2), N_int, j, & +! doubles, n_doubles) +! +! ! Compute Hij for all doubles +! ! --------------------------- +! +! m_b = psi_bilinear_matrix_order(l_a)+1 +! mcol = psi_bilinear_matrix_transp_columns(m_b) +! do j=1,n_doubles +! tmp_det2(1:N_int,2) = psi_det_beta_unique(1:N_int, doubles(j) ) +! call i_H_j_double_alpha_beta(tmp_det,tmp_det2,N_int,hij) +! do while (mcol /= doubles(j)) +! m_b = m_b+1 +! if (m_b > N_det) exit +! mcol = psi_bilinear_matrix_transp_columns(m_b) +! enddo +! m_a = psi_bilinear_matrix_order_reverse(m_b) ! v_0(m_a, 1:N_st) += hij * psi_bilinear_matrix_values(k_a,1:N_st) ! v_0(k_a, 1:N_st) += hij * psi_bilinear_matrix_values(m_a,1:N_st) - enddo - - endif - l_a = l_a+1 - if (l_a > N_det) exit - lrow = psi_bilinear_matrix_rows(l_a) - lcol = psi_bilinear_matrix_columns(l_a) - enddo - - k_b = k_b+1 - if (k_b > N_det) exit - kcol = psi_bilinear_matrix_transp_columns(k_b) - enddo - - do j=1,n_singles - is_single_a( singles(j) ) = .False. - enddo - - enddo +! enddo +! +! endif +! l_a = l_a+1 +! if (l_a > N_det) exit +! lrow = psi_bilinear_matrix_rows(l_a) +! lcol = psi_bilinear_matrix_columns(l_a) +! enddo +! +! k_b = k_b+1 +! if (k_b > N_det) exit +! kcol = psi_bilinear_matrix_transp_columns(k_b) +! enddo +! +! do j=1,n_singles +! is_single_a( singles(j) ) = .False. +! enddo +! +! enddo end @@ -946,7 +1027,7 @@ subroutine H_S2_u_0_nstates_test(v_0,s_0,u_0,H_jj,S2_jj,n,keys_tmp,Nint,N_st,sze integer :: degree integer :: exc(0:2,2,2) call get_excitation(keys_tmp(1,1,j),keys_tmp(1,1,i),exc,degree,phase,Nint) - if ((degree == 2).and.(exc(0,1,1)==1)) cycle +! if ((degree == 2).and.(exc(0,1,1)==1)) cycle ! if ((degree > 1)) cycle ! if (exc(0,1,2) /= 0) cycle call i_H_j(keys_tmp(1,1,j),keys_tmp(1,1,i),Nint,hij) diff --git a/src/Determinants/slater_rules.irp.f b/src/Determinants/slater_rules.irp.f index 56ad5617..1e0cb0a8 100644 --- a/src/Determinants/slater_rules.irp.f +++ b/src/Determinants/slater_rules.irp.f @@ -2566,13 +2566,14 @@ subroutine i_H_j_double_alpha_beta(key_i,key_j,Nint,hij) double precision, intent(out) :: hij integer :: exc(0:2,2,2) - double precision :: phase + double precision :: phase, phase2 double precision, external :: get_mo_bielec_integral PROVIDE big_array_exchange_integrals mo_bielec_integrals_in_map call get_mono_excitation_spin(key_i(1,1),key_j(1,1),exc(0,1,1),phase,Nint) - call get_mono_excitation_spin(key_i(1,2),key_j(1,2),exc(0,1,2),phase,Nint) + call get_mono_excitation_spin(key_i(1,2),key_j(1,2),exc(0,1,2),phase2,Nint) + phase = phase*phase2 if (exc(1,1,1) == exc(1,2,2)) then 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 diff --git a/src/Determinants/spindeterminants.irp.f b/src/Determinants/spindeterminants.irp.f index 4bb35979..783474f9 100644 --- a/src/Determinants/spindeterminants.irp.f +++ b/src/Determinants/spindeterminants.irp.f @@ -389,6 +389,7 @@ BEGIN_PROVIDER [ double precision, psi_bilinear_matrix_values, (N_det,N_states) &BEGIN_PROVIDER [ integer, psi_bilinear_matrix_rows , (N_det) ] &BEGIN_PROVIDER [ integer, psi_bilinear_matrix_columns, (N_det) ] &BEGIN_PROVIDER [ integer, psi_bilinear_matrix_order , (N_det) ] +&BEGIN_PROVIDER [ integer, psi_bilinear_matrix_columns_loc, (N_det_beta_unique+1) ] use bitmasks implicit none BEGIN_DOC @@ -428,6 +429,17 @@ BEGIN_PROVIDER [ double precision, psi_bilinear_matrix_values, (N_det,N_states) do l=1,N_states call dset_order(psi_bilinear_matrix_values(1,l),psi_bilinear_matrix_order,N_det) enddo + psi_bilinear_matrix_columns_loc(1:N_det_beta_unique) = -1 + psi_bilinear_matrix_columns_loc(1) = 1 + do k=2,N_det + if (psi_bilinear_matrix_columns(k) == psi_bilinear_matrix_columns(k-1)) then + cycle + else + l = psi_bilinear_matrix_columns(k) + psi_bilinear_matrix_columns_loc(l) = k + endif + enddo + psi_bilinear_matrix_columns_loc(N_det_beta_unique+1) = N_det+1 deallocate(to_sort) END_PROVIDER