LCPQ/quantum_package
10
0
Fork 0
mirror of https://github.com/LCPQ/quantum_package synced 2024-12-22 20:35:19 +01:00
Code Issues Releases Wiki Activity

Removed spin contaminants of Davidson

Browse Source
This commit is contained in:
Anthony Scemama 2016-11-18 15:06:33 +01:00
parent 278c961c0f
commit 5e3201cea9
2 changed files with 59 additions and 24 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

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

@ -781,27 +781,40 @@ 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 if (s2_eig) then
logical :: state_ok(N_st_diag*davidson_sze_max) logical :: state_ok(N_st_diag*davidson_sze_max)
do k=1,shift2 do k=1,shift2
state_ok(k) = (dabs(s2(k)-expected_s2) < 1.d0) state_ok(k) = (dabs(s2(k)-expected_s2) < 0.5d0)
enddo enddo
do k=1,shift2 else
if (.not. state_ok(k)) then state_ok(k) = .True.
do l=k+1,shift2
if (state_ok(l)) then
call dswap(shift2, y(1,k), 1, y(1,l), 1)
call dswap(1, s2(k), 1, s2(l), 1)
call dswap(1, lambda(k), 1, lambda(l), 1)
state_ok(k) = .True.
state_ok(l) = .False.
exit
endif
enddo
endif
enddo
endif endif
do k=1,shift2
if (.not. state_ok(k)) then
do l=k+1,shift2
if (state_ok(l)) then
call dswap(shift2, y(1,k), 1, y(1,l), 1)
call dswap(1, s2(k), 1, s2(l), 1)
call dswap(1, lambda(k), 1, lambda(l), 1)
state_ok(k) = .True.
state_ok(l) = .False.
exit
endif
enddo
endif
! Randomize components with bad <S2>
if (.not. state_ok(k)) then
do i=1,shift2
call random_number(r1)
call random_number(r2)
r1 = dsqrt(-2.d0*dlog(r1))
r2 = dtwo_pi*r2
y(i,k) = r1*dcos(r2)
lambda(k) = 1.d0
enddo
endif
enddo
! ! Compute overlap with U_in ! ! Compute overlap with U_in
! ! ------------------------- ! ! -------------------------
@ -852,11 +865,22 @@ subroutine davidson_diag_hjj_sjj_mrcc(dets_in,u_in,H_jj,S2_jj,energies,dim_in,sz
! ----------------------- ! -----------------------
do k=1,N_st_diag do k=1,N_st_diag
do i=1,sze do i=1,sze
U(i,shift2+k) = (lambda(k) * U(i,shift2+k) - W(i,shift2+k) ) & 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 & * (1.d0 + s2(k) * U(i,shift2+k) - S(i,shift2+k) - S_z2_Sz &
)/max(H_jj(i) - lambda (k),1.d-2) )/max(H_jj(i) - lambda (k),1.d-2)
enddo enddo
! else
! ! Randomize components with bad <S2>
! do i=1,sze
! call random_number(r1)
! call random_number(r2)
! r1 = dsqrt(-2.d0*dlog(r1))
! r2 = dtwo_pi*r2
! U(i,shift2+k) = r1*dcos(r2)
! enddo
! endif
if (k <= N_st) then if (k <= N_st) then
residual_norm(k) = u_dot_u(U(1,shift2+k),sze) residual_norm(k) = u_dot_u(U(1,shift2+k),sze)
to_print(1,k) = lambda(k) + nuclear_repulsion to_print(1,k) = lambda(k) + nuclear_repulsion

11
src/Davidson/diagonalization_hs2.irp.f
View File

@ -587,6 +587,17 @@ subroutine davidson_diag_hjj_sjj_mmap(dets_in,u_in,H_jj,S2_jj,energies,dim_in,sz
enddo enddo
endif endif
enddo enddo
! Randomize components with bad <S2>
if (.not. state_ok(k)) then
do i=1,shift2
call random_number(r1)
call random_number(r2)
r1 = dsqrt(-2.d0*dlog(r1))
r2 = dtwo_pi*r2
y(i,k) = r1*dcos(r2)
lambda(k) = 1.d0
enddo
endif
endif endif