From 0f2bdedb902c14c64defa6e5e50577332599b4f3 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Thu, 9 Apr 2015 22:32:25 +0200 Subject: [PATCH] Moves psi_cas in Dets --- src/MRCC/mrcc.irp.f | 27 +++---- src/MRCC/mrcc_dress.irp.f | 89 +++++++++++++---------- src/MRCC/mrcc_utils.irp.f | 145 +++----------------------------------- 3 files changed, 70 insertions(+), 191 deletions(-) diff --git a/src/MRCC/mrcc.irp.f b/src/MRCC/mrcc.irp.f index 48587b6d..71b87ddb 100644 --- a/src/MRCC/mrcc.irp.f +++ b/src/MRCC/mrcc.irp.f @@ -10,10 +10,6 @@ end subroutine run implicit none - print *, N_det - print *, N_det_cas - print *, N_det_sd - integer :: i,j print *, 'CAS' print *, '===' @@ -22,13 +18,13 @@ subroutine run call debug_det(psi_cas(1,1,i),N_int) enddo - print *, 'SD' - print *, '==' - do i=1,N_det_sd - print *, psi_sd_coefs(i,:) - call debug_det(psi_sd(1,1,i),N_int) - enddo - print *, 'xxx', 'Energy CAS+SD', ci_energy +! print *, 'SD' +! print *, '==' +! do i=1,N_det_non_cas +! print *, psi_non_cas_coefs(i,:) +! call debug_det(psi_non_cas(1,1,i),N_int) +! enddo + call write_double(6,ci_energy(1),"Initial CI energy") end subroutine run_mrcc_test implicit none @@ -46,14 +42,7 @@ end subroutine run_mrcc implicit none integer :: i,j - print *, 'MRCC' - print *, '====' - print *, '' - print *, 'CAS+SD energy : ', ci_energy_dressed(:) - print *, '' -! call diagonalize_ci_dressed -! call save_wavefunction_unsorted double precision :: E_new, E_old, delta_e integer :: iteration E_new = 0.d0 @@ -71,5 +60,7 @@ subroutine run_mrcc delta_E = dabs(E_new - E_old) call write_double(6,ci_energy_dressed(1),"MRCC energy") enddo + call ezfio_set_mrcc_energy(ci_energy_dressed(1)) +! call save_wavefunction end diff --git a/src/MRCC/mrcc_dress.irp.f b/src/MRCC/mrcc_dress.irp.f index 2e913b02..96f278ca 100644 --- a/src/MRCC/mrcc_dress.irp.f +++ b/src/MRCC/mrcc_dress.irp.f @@ -28,38 +28,44 @@ subroutine mrcc_dress(delta_ij_,Ndet,i_generator,n_selected,det_buffer,Nint,ipro integer :: h1,h2,p1,p2,s1,s2 integer(bit_kind):: tmp_det(Nint,2) integer :: iint, ipos -! integer :: istate, i_sd, i_cas + integer :: i_state, k_sd, l_sd, i_I, i_alpha ! |I> ! |alpha> - do i=1,N_tq - call get_excitation_degree_vector(psi_sd,tq(1,1,i),degree_alpha,Nint,N_det_sd,idx_alpha) + do i_alpha=1,N_tq + call get_excitation_degree_vector(psi_non_cas,tq(1,1,i_alpha),degree_alpha,Nint,N_det_non_cas,idx_alpha) ! |I> - do j=1,N_det_cas + do i_I=1,N_det_cas ! Find triples and quadruple grand parents - call get_excitation_degree(tq(1,1,i),psi_cas(1,1,j),degree,Nint) + call get_excitation_degree(tq(1,1,i_alpha),psi_cas(1,1,i_I),degree,Nint) if (degree > 4) then cycle endif - dIa(:) = 0.d0 + + do i_state=1,N_states + dIa(i_state) = 0.d0 + enddo + ! |alpha> - do k=1,idx_alpha(0) - call get_excitation_degree(psi_cas(1,1,j),psi_sd(1,1,idx_alpha(k)),degree,Nint) + do k_sd=1,idx_alpha(0) + call get_excitation_degree(psi_cas(1,1,i_I),psi_non_cas(1,1,idx_alpha(k_sd)),degree,Nint) if (degree > 2) then cycle endif - ! + ! ! - call i_h_j(psi_cas(1,1,j),psi_sd(1,1,idx_alpha(k)),Nint,hIk) - dIk(:) = hIk * lambda_mrcc(idx_alpha(k),:) - ! Exc(k -> alpha) - call get_excitation(psi_sd(1,1,idx_alpha(k)),tq(1,1,i),exc,degree,phase,Nint) + call i_h_j(psi_cas(1,1,i_I),psi_non_cas(1,1,idx_alpha(k_sd)),Nint,hIk) + do i_state=1,N_states + dIk(i_state) = hIk * lambda_mrcc(idx_alpha(k_sd),i_state) + enddo + ! |l> = Exc(k -> alpha) |I> + call get_excitation(psi_non_cas(1,1,idx_alpha(k_sd)),tq(1,1,i_alpha),exc,degree,phase,Nint) call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2) - tmp_det(1:Nint,1:2) = psi_cas(1,1,j) + tmp_det(1:Nint,1:2) = psi_cas(1,1,i_I) ! Hole (see list_to_bitstring) iint = ishft(h1-1,-bit_kind_shift) + 1 ipos = h1-ishft((iint-1),bit_kind_shift)-1 @@ -81,28 +87,37 @@ subroutine mrcc_dress(delta_ij_,Ndet,i_generator,n_selected,det_buffer,Nint,ipro tmp_det(iint,s2) = ibset(tmp_det(iint,s2),ipos) endif - dka(:) = 0.d0 - do l=k+1,idx_alpha(0) - call get_excitation_degree(tmp_det,psi_sd(1,1,idx_alpha(l)),degree,Nint) + ! + do i_state=1,N_states + dka(i_state) = 0.d0 + enddo + do l_sd=k_sd+1,idx_alpha(0) + call get_excitation_degree(tmp_det,psi_non_cas(1,1,idx_alpha(l_sd)),degree,Nint) if (degree == 0) then - call get_excitation(psi_cas(1,1,j),psi_sd(1,1,idx_alpha(l)),exc,degree,phase2,Nint) - call i_h_j(psi_cas(1,1,j),psi_sd(1,1,idx_alpha(l)),Nint,hIl) - dka(:) = hIl * lambda_mrcc(idx_alpha(l),:) * phase * phase2 + call get_excitation(psi_cas(1,1,i_I),psi_non_cas(1,1,idx_alpha(l_sd)),exc,degree,phase2,Nint) + call i_h_j(psi_cas(1,1,i_I),psi_non_cas(1,1,idx_alpha(l_sd)),Nint,hIl) + do i_state=1,N_states + dka(i_state) = hIl * lambda_mrcc(idx_alpha(l_sd),i_state) * phase * phase2 + enddo exit endif enddo - do l=1,N_states - dIa(l) += dka(l)*dIk(l) + do i_state=1,N_states + dIa(i_state) = dIa(i_state) + dIk(i_state) * dka(i_state) enddo enddo - ci_inv(1:N_states) = 1.d0/psi_cas_coefs(j,1:N_states) - do l=1,idx_alpha(0) - k = idx_alpha(l) - call i_h_j(tq(1,1,i),psi_sd(1,1,idx_alpha(l)),Nint,hla) + + do i_state=1,N_states + ci_inv(i_state) = 1.d0/psi_cas_coefs(i_I,i_state) + enddo + + do l_sd=1,idx_alpha(0) + k_sd = idx_alpha(l_sd) + call i_h_j(tq(1,1,i_alpha),psi_non_cas(1,1,idx_alpha(l_sd)),Nint,hla) do m=1,N_states - delta_ij_(idx_sd(k),idx_cas(j),m) += dIa(m) * hla - delta_ij_(idx_cas(j),idx_sd(k),m) += dIa(m) * hla - delta_ij_(idx_cas(j),idx_cas(j),m) -= dIa(m) * hla * ci_inv(m) * psi_sd_coefs(k,m) + delta_ij_(idx_non_cas(k_sd),idx_cas(i_I),m) += dIa(m) * hla + delta_ij_(idx_cas(i_I),idx_non_cas(k_sd),m) += dIa(m) * hla + delta_ij_(idx_cas(i_I),idx_cas(i_I),m) -= dIa(m) * hla * ci_inv(m) * psi_non_cas_coefs(k_sd,m) enddo enddo enddo @@ -115,13 +130,13 @@ end -subroutine mrcc_dress_simple(delta_ij_sd_,Ndet_sd,i_generator,n_selected,det_buffer,Nint,iproc) +subroutine mrcc_dress_simple(delta_ij_non_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_sd - double precision, intent(inout) :: delta_ij_sd_(Ndet_sd,Ndet_sd,*) + integer, intent(in) :: Ndet_non_cas + double precision, intent(inout) :: delta_ij_non_cas_(Ndet_non_cas,Ndet_non_cas,*) integer(bit_kind), intent(in) :: det_buffer(Nint,2,n_selected) integer :: i,j,k,m @@ -143,18 +158,18 @@ subroutine mrcc_dress_simple(delta_ij_sd_,Ndet_sd,i_generator,n_selected,det_buf double precision :: f(N_states) do i=1,N_tq - call get_excitation_degree_vector(psi_sd,tq(1,1,i),degree,Nint,Ndet_sd,idx) + call get_excitation_degree_vector(psi_non_cas,tq(1,1,i),degree,Nint,Ndet_non_cas,idx) call i_h_j(tq(1,1,i),tq(1,1,i),Nint,haa) do m=1,N_states f(m) = 1.d0/(ci_electronic_energy(m)-haa) enddo do k=1,idx(0) - call i_h_j(tq(1,1,i),psi_sd(1,1,idx(k)),Nint,hka) + call i_h_j(tq(1,1,i),psi_non_cas(1,1,idx(k)),Nint,hka) do j=k,idx(0) - call i_h_j(tq(1,1,i),psi_sd(1,1,idx(j)),Nint,haj) + call i_h_j(tq(1,1,i),psi_non_cas(1,1,idx(j)),Nint,haj) do m=1,N_states - delta_ij_sd_(idx(k), idx(j),m) += haj*hka* f(m) - delta_ij_sd_(idx(j), idx(k),m) += haj*hka* f(m) + delta_ij_non_cas_(idx(k), idx(j),m) += haj*hka* f(m) + delta_ij_non_cas_(idx(j), idx(k),m) += haj*hka* f(m) enddo enddo enddo diff --git a/src/MRCC/mrcc_utils.irp.f b/src/MRCC/mrcc_utils.irp.f index a2b91645..4a1fc1cc 100644 --- a/src/MRCC/mrcc_utils.irp.f +++ b/src/MRCC/mrcc_utils.irp.f @@ -1,130 +1,3 @@ - use bitmasks - -BEGIN_PROVIDER [ integer(bit_kind), cas_bitmask, (N_int,2,N_cas_bitmask) ] - implicit none - BEGIN_DOC - ! Bitmasks for CAS reference determinants. (N_int, alpha/beta, CAS reference) - END_DOC - logical :: exists - integer :: i - PROVIDE ezfio_filename - - call ezfio_has_bitmasks_cas(exists) - if (exists) then - call ezfio_get_bitmasks_cas(cas_bitmask) - else - do i=1,N_cas_bitmask - cas_bitmask(:,:,i) = iand(not(HF_bitmask(:,:)),full_ijkl_bitmask(:,:)) - enddo - endif - -END_PROVIDER - -BEGIN_PROVIDER [ integer, N_det_cas ] - implicit none - BEGIN_DOC - ! Number of generator detetrminants - END_DOC - integer :: i,k,l - logical :: good - call write_time(output_dets) - N_det_cas = 0 - do i=1,N_det - do l=1,n_cas_bitmask - good = .True. - do k=1,N_int - good = good .and. ( & - iand(not(cas_bitmask(k,1,l)), psi_det(k,1,i)) == & - iand(not(cas_bitmask(k,1,l)), psi_det(k,1,1)) ) .and. ( & - iand(not(cas_bitmask(k,2,l)), psi_det(k,2,i)) == & - iand(not(cas_bitmask(k,2,l)), psi_det(k,2,1)) ) - enddo - if (good) then - exit - endif - enddo - if (good) then - N_det_cas += 1 - endif - enddo - N_det_cas = max(N_det_cas, 1) - call write_int(output_dets,N_det_cas, 'Number of determinants in the CAS') -END_PROVIDER - - BEGIN_PROVIDER [ integer(bit_kind), psi_cas, (N_int,2,N_det_cas) ] -&BEGIN_PROVIDER [ double precision, psi_cas_coefs, (N_det_cas,n_states) ] -&BEGIN_PROVIDER [ integer, idx_cas, (N_det_cas) ] - implicit none - BEGIN_DOC - ! For Single reference wave functions, the generator is the - ! Hartree-Fock determinant - END_DOC - integer :: i, k, l, m - logical :: good - m=0 - do i=1,N_det - do l=1,n_cas_bitmask - good = .True. - do k=1,N_int - good = good .and. ( & - iand(not(cas_bitmask(k,1,l)), psi_det(k,1,i)) == & - iand(not(cas_bitmask(k,1,l)), psi_det(k,1,1)) ) .and. ( & - iand(not(cas_bitmask(k,2,l)), psi_det(k,2,i)) == & - iand(not(cas_bitmask(k,2,l)), psi_det(k,2,1)) ) - enddo - if (good) then - exit - endif - enddo - if (good) then - m = m+1 - do k=1,N_int - psi_cas(k,1,m) = psi_det(k,1,i) - psi_cas(k,2,m) = psi_det(k,2,i) - enddo - idx_cas(m) = i - do k=1,N_states - psi_cas_coefs(m,k) = psi_coef(i,k) - enddo - endif - enddo - -END_PROVIDER - - - - - BEGIN_PROVIDER [ integer(bit_kind), psi_sd, (N_int,2,N_det) ] -&BEGIN_PROVIDER [ double precision, psi_sd_coefs, (N_det,n_states) ] -&BEGIN_PROVIDER [ integer, idx_sd, (N_det) ] -&BEGIN_PROVIDER [ integer, N_det_sd] - implicit none - BEGIN_DOC - ! SD - END_DOC - integer :: i_sd,j,k - integer :: degree - logical :: in_cas - i_sd =0 - do k=1,N_det - in_cas = .False. - do j=1,N_det_cas - call get_excitation_degree(psi_cas(1,1,j), psi_det(1,1,k), degree, N_int) - if (degree == 0) then - in_cas = .True. - exit - endif - enddo - if (.not.in_cas) then - double precision :: hij - i_sd += 1 - psi_sd(1:N_int,1:2,i_sd) = psi_det(1:N_int,1:2,k) - psi_sd_coefs(i_sd,1:N_states) = psi_coef(k,1:N_states) - idx_sd(i_sd) = k - endif - enddo - N_det_sd = i_sd -END_PROVIDER BEGIN_PROVIDER [ double precision, lambda_mrcc, (psi_det_size,n_states) ] implicit none @@ -133,15 +6,15 @@ BEGIN_PROVIDER [ double precision, lambda_mrcc, (psi_det_size,n_states) ] END_DOC integer :: i,k double precision :: ihpsi(N_states) - do i=1,N_det_sd - call i_h_psi(psi_sd(1,1,i), psi_cas, psi_cas_coefs, N_int, N_det_cas, & + do i=1,N_det_non_cas + call i_h_psi(psi_non_cas(1,1,i), psi_cas, psi_cas_coefs, N_int, N_det_cas, & size(psi_cas_coefs,1), n_states, ihpsi) double precision :: hij do k=1,N_states if (dabs(ihpsi(k)) < 1.d-6) then lambda_mrcc(i,k) = 0.d0 else - lambda_mrcc(i,k) = psi_sd_coefs(i,k)/ihpsi(k) + lambda_mrcc(i,k) = psi_non_cas_coefs(i,k)/ihpsi(k) lambda_mrcc(i,k) = min( lambda_mrcc (i,k),0.d0 ) endif enddo @@ -158,13 +31,13 @@ BEGIN_PROVIDER [ character*(32), dressing_type ] dressing_type = "MRCC" END_PROVIDER -BEGIN_PROVIDER [ double precision, delta_ij_sd, (N_det_sd, N_det_sd,N_states) ] +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_sd = 0.d0 - call H_apply_mrcc_simple(delta_ij_sd,N_det_sd) + 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) ] @@ -178,9 +51,9 @@ BEGIN_PROVIDER [ double precision, delta_ij, (N_det,N_det,N_states) ] call H_apply_mrcc(delta_ij,N_det) else if (dressing_type == "Simple") then do m=1,N_states - do j=1,N_det_sd - do i=1,N_det_sd - delta_ij(idx_sd(i),idx_sd(j),m) = delta_ij_sd(i,j,m) + 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