LCPQ/quantum_package
10
0
Fork 0
mirror of https://github.com/LCPQ/quantum_package synced 2024-06-29 16:34:50 +02:00
Code Issues Releases Wiki Activity

State following in MRCC

Browse Source
This commit is contained in:
Anthony Scemama 2016-11-15 11:23:00 +01:00
parent c366c201eb
commit ee658adeb7
2 changed files with 50 additions and 36 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

72
plugins/MRCC_Utils/davidson.irp.f
View File

@ -628,7 +628,7 @@ subroutine davidson_diag_hjj_sjj_mrcc(dets_in,u_in,H_jj,S2_jj,energies,dim_in,sz
integer :: k_pairs, kl integer :: k_pairs, kl
integer :: iter2 integer :: iter2
double precision, allocatable :: W(:,:), U(:,:), S(:,:) double precision, allocatable :: W(:,:), U(:,:), S(:,:), overlap(:,:)
double precision, allocatable :: y(:,:), h(:,:), lambda(:), s2(:) double precision, allocatable :: y(:,:), h(:,:), lambda(:), s2(:)
double precision, allocatable :: c(:), s_(:,:), s_tmp(:,:) double precision, allocatable :: c(:), s_(:,:), s_tmp(:,:)
double precision :: diag_h_mat_elem double precision :: diag_h_mat_elem
@ -688,16 +688,17 @@ subroutine davidson_diag_hjj_sjj_mrcc(dets_in,u_in,H_jj,S2_jj,energies,dim_in,sz
itermax = min(davidson_sze_max, sze/N_st_diag) itermax = min(davidson_sze_max, sze/N_st_diag)
allocate( & allocate( &
W(sze_8,N_st_diag*itermax), & W(sze_8,N_st_diag*itermax), &
U(sze_8,N_st_diag*itermax), & U(sze_8,N_st_diag*itermax), &
S(sze_8,N_st_diag*itermax), & S(sze_8,N_st_diag*itermax), &
h(N_st_diag*itermax,N_st_diag*itermax), & h(N_st_diag*itermax,N_st_diag*itermax), &
y(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_(N_st_diag*itermax,N_st_diag*itermax), &
s_tmp(N_st_diag*itermax,N_st_diag*itermax), & s_tmp(N_st_diag*itermax,N_st_diag*itermax), &
residual_norm(N_st_diag), & residual_norm(N_st_diag), &
c(N_st_diag*itermax), & c(N_st_diag*itermax), &
s2(N_st_diag*itermax), & s2(N_st_diag*itermax), &
overlap(N_st_diag*itermax,N_st_diag*itermax), &
lambda(N_st_diag*itermax)) lambda(N_st_diag*itermax))
h = 0.d0 h = 0.d0
@ -795,26 +796,39 @@ subroutine davidson_diag_hjj_sjj_mrcc(dets_in,u_in,H_jj,S2_jj,energies,dim_in,sz
s2(k) = s_(k,k) + S_z2_Sz s2(k) = s_(k,k) + S_z2_Sz
enddo enddo
if (s2_eig) then ! Compute overlap with U_in
logical :: state_ok(N_st_diag*davidson_sze_max) ! -------------------------
do k=1,shift2
state_ok(k) = (dabs(s2(k)-expected_s2) < 0.6d0) integer :: coord(2), order(N_st)
overlap = -1.d0
do k=1,N_st
do i=1,shift2
overlap(i,k) = dabs(y(i,k))
enddo enddo
do k=1,shift2 enddo
if (.not. state_ok(k)) then do k=1,N_st
do l=k+1,shift2 coord = maxloc(overlap)
if (state_ok(l)) then order( coord(2) ) = coord(1)
call dswap(shift2, y(1,k), 1, y(1,l), 1) overlap(coord(1),coord(2)) = -1.d0
call dswap(1, s2(k), 1, s2(l), 1) enddo
call dswap(1, lambda(k), 1, lambda(l), 1) overlap = y
state_ok(k) = .True. do k=1,N_st
state_ok(l) = .False. l = order(k)
exit if (k /= l) then
endif y(1:shift2,k) = overlap(1:shift2,l)
enddo endif
endif enddo
enddo do k=1,N_st
endif overlap(k,1) = lambda(k)
overlap(k,2) = s2(k)
enddo
do k=1,N_st
l = order(k)
if (k /= l) then
lambda(k) = overlap(l,1)
s2(k) = overlap(l,2)
endif
enddo
! Express eigenvectors of h in the determinant basis ! Express eigenvectors of h in the determinant basis

14
plugins/MRCC_Utils/mrcc_utils.irp.f
View File

@ -149,7 +149,7 @@ END_PROVIDER
allocate (eigenvectors(size(CI_eigenvectors_dressed,1),size(CI_eigenvectors_dressed,2)), & allocate (eigenvectors(size(CI_eigenvectors_dressed,1),size(CI_eigenvectors_dressed,2)), &
eigenvalues(size(CI_electronic_energy_dressed,1))) eigenvalues(size(CI_electronic_energy_dressed,1)))
do mrcc_state=N_states,1,-1 do mrcc_state=1,N_states
do j=1,min(N_states,N_det) do j=1,min(N_states,N_det)
do i=1,N_det do i=1,N_det
eigenvectors(i,j) = psi_coef(i,j) eigenvectors(i,j) = psi_coef(i,j)
@ -161,10 +161,12 @@ END_PROVIDER
output_determinants,mrcc_state) output_determinants,mrcc_state)
CI_eigenvectors_dressed(1:N_det,mrcc_state) = eigenvectors(1:N_det,mrcc_state) CI_eigenvectors_dressed(1:N_det,mrcc_state) = eigenvectors(1:N_det,mrcc_state)
CI_electronic_energy_dressed(mrcc_state) = eigenvalues(mrcc_state) CI_electronic_energy_dressed(mrcc_state) = eigenvalues(mrcc_state)
enddo if (mrcc_state == 1) then
do mrcc_state=N_states+1,N_states_diag do mrcc_state=N_states+1,N_states_diag
CI_eigenvectors_dressed(1:N_det,mrcc_state) = eigenvectors(1:N_det,mrcc_state) CI_eigenvectors_dressed(1:N_det,mrcc_state) = eigenvectors(1:N_det,mrcc_state)
CI_electronic_energy_dressed(mrcc_state) = eigenvalues(mrcc_state) CI_electronic_energy_dressed(mrcc_state) = eigenvalues(mrcc_state)
enddo
endif
enddo enddo
call u_0_S2_u_0(CI_eigenvectors_s2_dressed,CI_eigenvectors_dressed,N_det,psi_det,N_int,& 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)) N_states_diag,size(CI_eigenvectors_dressed,1))
@ -685,7 +687,6 @@ END_PROVIDER
if(is_active_exc(pp)) cycle if(is_active_exc(pp)) cycle
lref = 0 lref = 0
AtB(pp) = 0.d0 AtB(pp) = 0.d0
X(pp) = 0.d0
do II=1,N_det_ref do II=1,N_det_ref
call apply_hole_local(psi_ref(1,1,II), hh_exists(1, hh), myMask, ok, N_int) call apply_hole_local(psi_ref(1,1,II), hh_exists(1, hh), myMask, ok, N_int)
if(.not. ok) cycle if(.not. ok) cycle
@ -695,7 +696,6 @@ END_PROVIDER
if(ind == -1) cycle if(ind == -1) cycle
ind = psi_non_ref_sorted_idx(ind) ind = psi_non_ref_sorted_idx(ind)
call get_phase(myDet(1,1), psi_ref(1,1,II), phase, N_int) call get_phase(myDet(1,1), psi_ref(1,1,II), phase, N_int)
X(pp) = X(pp) + psi_ref_coef(II,s)*psi_ref_coef(II,s)
AtB(pp) += psi_non_ref_coef(ind, s) * psi_ref_coef(II, s) * phase AtB(pp) += psi_non_ref_coef(ind, s) * psi_ref_coef(II, s) * phase
lref(II) = ind lref(II) = ind
if(phase < 0.d0) lref(II) = -ind if(phase < 0.d0) lref(II) = -ind