From fb1e97ce1a620470a1032721b676676f038c8c56 Mon Sep 17 00:00:00 2001 From: Manu Date: Fri, 3 Jul 2015 18:50:07 +0200 Subject: [PATCH] bug in mrcc with tiny matrix for the dressing --- plugins/MRCC/H_apply.irp.f | 26 +++++------ plugins/MRCC/davidson.irp.f | 39 +++++++++++------ plugins/MRCC/mrcc_dress.irp.f | 16 +++---- plugins/MRCC/mrcc_utils.irp.f | 81 ++++++++++++++--------------------- 4 files changed, 77 insertions(+), 85 deletions(-) diff --git a/plugins/MRCC/H_apply.irp.f b/plugins/MRCC/H_apply.irp.f index cadaca32..7e479529 100644 --- a/plugins/MRCC/H_apply.irp.f +++ b/plugins/MRCC/H_apply.irp.f @@ -2,18 +2,20 @@ use bitmasks BEGIN_SHELL [ /usr/bin/env python ] from generate_h_apply import * -s = H_apply("mrcc_simple") -s.data["parameters"] = ", delta_ij_sd_, Ndet_sd" +s = H_apply("mrcc") +s.data["parameters"] = ", delta_ij_, delta_ii_,Ndet_cas, Ndet_non_cas" s.data["declarations"] += """ - integer, intent(in) :: Ndet_sd - double precision, intent(in) :: delta_ij_sd_(Ndet_sd,Ndet_sd,*) + integer, intent(in) :: Ndet_cas,Ndet_non_cas + double precision, intent(in) :: delta_ij_(Ndet_cas,Ndet_non_cas,*) + double precision, intent(in) :: delta_ii_(Ndet_cas,*) """ -s.data["keys_work"] = "call mrcc_dress_simple(delta_ij_sd_,Ndet_sd,i_generator,key_idx,keys_out,N_int,iproc)" -s.data["params_post"] += ", delta_ij_sd_, Ndet_sd" -s.data["params_main"] += "delta_ij_sd_, Ndet_sd" +s.data["keys_work"] = "call mrcc_dress(delta_ij_,delta_ii_,Ndet_cas,Ndet_non_cas,i_generator,key_idx,keys_out,N_int,iproc)" +s.data["params_post"] += ", delta_ij_, delta_ii_, Ndet_cas, Ndet_non_cas" +s.data["params_main"] += "delta_ij_, delta_ii_, Ndet_cas, Ndet_non_cas" s.data["decls_main"] += """ - integer, intent(in) :: Ndet_sd - double precision, intent(in) :: delta_ij_sd_(Ndet_sd,Ndet_sd,*) + integer, intent(in) :: Ndet_cas,Ndet_non_cas + double precision, intent(in) :: delta_ij_(Ndet_cas,Ndet_non_cas,*) + double precision, intent(in) :: delta_ii_(Ndet_cas,*) """ s.data["finalization"] = "" s.data["copy_buffer"] = "" @@ -22,11 +24,5 @@ s.data["size_max"] = "3072" print s - - -s.data["subroutine"] = "H_apply_mrcc" -s.data["keys_work"] = "call mrcc_dress(delta_ij_sd_,Ndet_sd,i_generator,key_idx,keys_out,N_int,iproc)" -print s - END_SHELL diff --git a/plugins/MRCC/davidson.irp.f b/plugins/MRCC/davidson.irp.f index 57900082..d30591a3 100644 --- a/plugins/MRCC/davidson.irp.f +++ b/plugins/MRCC/davidson.irp.f @@ -35,11 +35,16 @@ subroutine davidson_diag_mrcc(dets_in,u_in,energies,dim_in,sze,N_st,Nint,iunit,i allocate(H_jj(sze)) !$OMP PARALLEL DEFAULT(NONE) & - !$OMP SHARED(sze,H_jj,dets_in,Nint,istate,delta_ij) & + !$OMP SHARED(sze,H_jj,N_det_cas,dets_in,Nint,istate,delta_ij,delta_ii,idx_cas) & !$OMP PRIVATE(i) !$OMP DO SCHEDULE(guided) do i=1,sze - H_jj(i) = diag_h_mat_elem(dets_in(1,1,i),Nint) + delta_ij(i,i,istate) + 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_cas + H_jj(idx_cas(i)) += delta_ii(i,istate) enddo !$OMP END DO !$OMP END PARALLEL @@ -370,16 +375,16 @@ subroutine H_u_0_mrcc(v_0,u_0,H_jj,n,keys_tmp,Nint,istate) integer, allocatable :: idx(:) double precision :: hij double precision, allocatable :: vt(:) - integer :: i,j,k,l, jj + integer :: i,j,k,l, jj,ii integer :: i0, j0 ASSERT (Nint > 0) ASSERT (Nint == N_int) ASSERT (n>0) - PROVIDE ref_bitmask_energy delta_ij + PROVIDE ref_bitmask_energy delta_ij delta_ii integer, parameter :: block_size = 157 !$OMP PARALLEL DEFAULT(NONE) & - !$OMP PRIVATE(i,hij,j,k,idx,jj,vt) & - !$OMP SHARED(n,H_jj,u_0,keys_tmp,Nint,v_0,istate,delta_ij) + !$OMP PRIVATE(i,hij,j,k,idx,jj,ii,vt) & + !$OMP SHARED(n_det_cas,n_det_non_cas,idx_cas,idx_non_cas,n,H_jj,u_0,keys_tmp,Nint,v_0,istate,delta_ij) !$OMP DO SCHEDULE(static) do i=1,n v_0(i) = H_jj(i) * u_0(i) @@ -389,20 +394,30 @@ subroutine H_u_0_mrcc(v_0,u_0,H_jj,n,keys_tmp,Nint,istate) Vt = 0.d0 !$OMP DO SCHEDULE(guided) do i=1,n -! idx(0) = i -! call filter_connected_davidson(keys_tmp,keys_tmp(1,1,i),Nint,i-1,idx) -! do jj=1,idx(0) -! j = idx(jj) + idx(0) = i + call filter_connected_davidson(keys_tmp,keys_tmp(1,1,i),Nint,i-1,idx) + do jj=1,idx(0) + j = idx(jj) ! if ( (dabs(u_0(j)) > 1.d-7).or.((dabs(u_0(i)) > 1.d-7)) ) then - do j = 1, i-1 call i_H_j(keys_tmp(1,1,j),keys_tmp(1,1,i),Nint,hij) - hij = hij + delta_ij(j,i,istate) + hij = hij vt (i) = vt (i) + hij*u_0(j) vt (j) = vt (j) + hij*u_0(i) ! endif enddo enddo !$OMP END DO + + !$OMP DO SCHEDULE(guided) + do ii=1,n_det_cas + i = idx_cas(ii) + do jj = 1, n_det_non_cas + j = idx_non_cas(jj) + vt (i) = vt (i) + delta_ij(j,i,istate)*u_0(j) + vt (j) = vt (j) + delta_ij(j,i,istate)*u_0(i) + enddo + enddo + !$OMP END DO !$OMP CRITICAL do i=1,n v_0(i) = v_0(i) + vt(i) diff --git a/plugins/MRCC/mrcc_dress.irp.f b/plugins/MRCC/mrcc_dress.irp.f index 23c15d13..53f94f5b 100644 --- a/plugins/MRCC/mrcc_dress.irp.f +++ b/plugins/MRCC/mrcc_dress.irp.f @@ -12,13 +12,14 @@ BEGIN_PROVIDER [ integer(omp_lock_kind), psi_cas_lock, (psi_det_size) ] END_PROVIDER -subroutine mrcc_dress(delta_ij_,Ndet,i_generator,n_selected,det_buffer,Nint,iproc) +subroutine mrcc_dress(delta_ij_, delta_ii_, Ndet_cas, Ndet_non_cas,i_generator,n_selected,det_buffer,Nint,iproc) use bitmasks implicit none integer, intent(in) :: i_generator,n_selected, Nint, iproc - integer, intent(in) :: Ndet - double precision, intent(inout) :: delta_ij_(Ndet,Ndet,*) + integer, intent(in) :: Ndet_cas, Ndet_non_cas + double precision, intent(inout) :: delta_ij_(Ndet_cas,Ndet_non_cas,*) + double precision, intent(inout) :: delta_ii_(Ndet_cas,*) integer(bit_kind), intent(in) :: det_buffer(Nint,2,n_selected) integer :: i,j,k,l @@ -42,7 +43,7 @@ subroutine mrcc_dress(delta_ij_,Ndet,i_generator,n_selected,det_buffer,Nint,ipro call find_triples_and_quadruples(i_generator,n_selected,det_buffer,Nint,tq,N_tq) - allocate (dIa_hla(N_states,Ndet)) + allocate (dIa_hla(N_states,Ndet_non_cas)) ! |I> @@ -136,12 +137,11 @@ subroutine mrcc_dress(delta_ij_,Ndet,i_generator,n_selected,det_buffer,Nint,ipro do l_sd=1,idx_alpha(0) k_sd = idx_alpha(l_sd) do i_state=1,N_states - delta_ij_(idx_non_cas(k_sd),idx_cas(i_I),i_state) += dIa_hla(i_state,k_sd) - delta_ij_(idx_cas(i_I),idx_non_cas(k_sd),i_state) += dIa_hla(i_state,k_sd) + delta_ij_(i_I,k_sd,i_state) += dIa_hla(i_state,k_sd) if(dabs(psi_cas_coef(i_I,i_state)).ge.5.d-5)then - delta_ij_(idx_cas(i_I),idx_cas(i_I),i_state) -= dIa_hla(i_state,k_sd) * ci_inv(i_state) * psi_non_cas_coef(k_sd,i_state) + delta_ii_(i_I,i_state) -= dIa_hla(i_state,k_sd) * ci_inv(i_state) * psi_non_cas_coef(k_sd,i_state) else - delta_ij_(idx_cas(i_I),idx_cas(i_I),i_state) = 0.d0 + delta_ii_(i_I,i_state) = 0.d0 endif enddo enddo diff --git a/plugins/MRCC/mrcc_utils.irp.f b/plugins/MRCC/mrcc_utils.irp.f index b175d655..016f3762 100644 --- a/plugins/MRCC/mrcc_utils.irp.f +++ b/plugins/MRCC/mrcc_utils.irp.f @@ -12,13 +12,11 @@ do i=1,N_det_non_cas size(psi_cas_coef,1), n_states, ihpsi) call i_h_j(psi_non_cas(1,1,i),psi_non_cas(1,1,i),N_int,hii) do k=1,N_states - lambda_pert(k,i) = 1.d0 / (psi_cas_energy_diagonalized(k)-hii) - lambda_mrcc(k,i) = psi_non_cas_coef(i,k)/ihpsi(k) if (dabs(ihpsi(k)).le.1.d-3) then - lambda_mrcc(k,i) = 1.d0 / (psi_cas_energy_diagonalized(k)-hii) - icount_manu = icount_manu+1 - cycle + lambda_mrcc(k,i) = lambda_pert(k,i) + else + lambda_mrcc(k,i) = psi_non_cas_coef(i,k)/ihpsi(k) endif enddo enddo @@ -28,66 +26,49 @@ END_PROVIDER +!BEGIN_PROVIDER [ double precision, delta_ij_non_cas, (N_det_non_cas, N_det_non_cas,N_states) ] +!implicit none +!BEGIN_DOC +!! Dressing matrix in SD basis +!END_DOC +!delta_ij_non_cas = 0.d0 +!call H_apply_mrcc_simple(delta_ij_non_cas,N_det_non_cas) +!END_PROVIDER -BEGIN_PROVIDER [ character*(32), dressing_type ] - implicit none - BEGIN_DOC - ! [ Simple | MRCC ] - END_DOC - dressing_type = "MRCC" -END_PROVIDER - -BEGIN_PROVIDER [ double precision, delta_ij_non_cas, (N_det_non_cas, N_det_non_cas,N_states) ] - implicit none - BEGIN_DOC - ! Dressing matrix in SD basis - END_DOC - delta_ij_non_cas = 0.d0 - call H_apply_mrcc_simple(delta_ij_non_cas,N_det_non_cas) -END_PROVIDER - -BEGIN_PROVIDER [ double precision, delta_ij, (N_det,N_det,N_states) ] + BEGIN_PROVIDER [ double precision, delta_ij, (N_det_cas,N_det_non_cas,N_states) ] +&BEGIN_PROVIDER [ double precision, delta_ii, (N_det_cas,N_states) ] implicit none BEGIN_DOC ! Dressing matrix in N_det basis END_DOC integer :: i,j,m delta_ij = 0.d0 - if (dressing_type == "MRCC") then - call H_apply_mrcc(delta_ij,N_det) - else if (dressing_type == "Simple") then - do m=1,N_states - do j=1,N_det_non_cas - do i=1,N_det_non_cas - delta_ij(idx_non_cas(i),idx_non_cas(j),m) = delta_ij_non_cas(i,j,m) - enddo - enddo - enddo - endif - do i = 1, N_det - do j = 1, N_det - do m = 1, N_states - if(isnan(delta_ij(j,i,m)))then - delta_ij(j,i,m) = 0.d0 - endif - enddo - enddo - enddo - + delta_ii = 0.d0 + call H_apply_mrcc(delta_ij,delta_ii,N_det_cas,N_det_non_cas) END_PROVIDER -BEGIN_PROVIDER [ double precision, h_matrix_dressed, (N_det,N_det) ] +BEGIN_PROVIDER [ double precision, h_matrix_dressed, (N_det,N_det,N_states) ] implicit none BEGIN_DOC ! Dressed H with Delta_ij END_DOC - integer :: i, j - do j=1,N_det - do i=1,N_det - h_matrix_dressed(i,j) = h_matrix_all_dets(i,j) + delta_ij(i,j,1) + integer :: i, j,istate,ii,jj + do istate = 1,N_states + do j=1,N_det + do i=1,N_det + h_matrix_dressed(i,j,istate) = h_matrix_all_dets(i,j) + enddo enddo + do ii = 1, N_det_cas + i =idx_cas(ii) + h_matrix_dressed(i,i,istate) += delta_ii(ii,istate) + do jj = 1, N_det_non_cas + j =idx_cas(jj) + h_matrix_dressed(i,j,istate) += delta_ij(ii,jj,istate) + h_matrix_dressed(j,i,istate) += delta_ij(ii,jj,istate) + enddo + enddo enddo - END_PROVIDER