LCPQ/qp2
9
1
Fork 0
mirror of https://github.com/QuantumPackage/qp2.git synced 2024-12-21 11:03:29 +01:00
Code Releases Activity

fixed conflict after TC S^2 merge

Browse Source
This commit is contained in:
AbdAmmar 2023-04-10 20:41:16 +02:00
parent d67861342a
commit 16955745f4
12 changed files with 1754 additions and 109 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
external/qp2-dependencies vendored

@ -1 +1 @@
Subproject commit ce14f57b50511825a9fedb096749200779d3f4d4 Subproject commit 6e23ebac001acae91d1c762ca934e09a9b7d614a

3
src/davidson/diagonalization_hs2_dressed.irp.f
View File

@ -461,7 +461,8 @@ subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,s2_out,energies,dim_in,sze,N_
integer :: lwork, info integer :: lwork, info
double precision, allocatable :: work(:) double precision, allocatable :: work(:)
y = h ! y = h
y = h_p
lwork = -1 lwork = -1
allocate(work(1)) allocate(work(1))
call dsygv(1,'V','U',shift2,y,size(y,1), & call dsygv(1,'V','U',shift2,y,size(y,1), &

103
src/determinants/h_apply.irp.f
View File

@ -69,9 +69,15 @@ subroutine resize_H_apply_buffer(new_size,iproc)
END_DOC END_DOC
PROVIDE H_apply_buffer_allocated PROVIDE H_apply_buffer_allocated
ASSERT (new_size > 0) ASSERT (new_size > 0)
ASSERT (iproc >= 0) ASSERT (iproc >= 0)
ASSERT (iproc < nproc) ASSERT (iproc < nproc)
if (N_det < 0) call abort() !irp_here//': N_det < 0')
if (N_int <= 0) call abort() !irp_here//': N_int <= 0')
if (new_size <= 0) call abort() !irp_here//': new_size <= 0')
if (iproc < 0) call abort() !irp_here//': iproc < 0')
if (iproc >= nproc) call abort() !irp_here//': iproc >= nproc')
allocate ( buffer_det(N_int,2,new_size), & allocate ( buffer_det(N_int,2,new_size), &
buffer_coef(new_size,N_states), & buffer_coef(new_size,N_states), &
@ -126,31 +132,34 @@ subroutine copy_H_apply_buffer_to_wf
ASSERT (N_int > 0) ASSERT (N_int > 0)
ASSERT (N_det > 0) ASSERT (N_det >= 0)
allocate ( buffer_det(N_int,2,N_det), buffer_coef(N_det,N_states) ) N_det_old = N_det
if (N_det > 0) then
allocate ( buffer_det(N_int,2,N_det), buffer_coef(N_det,N_states) )
! Backup determinants ! Backup determinants
j=0 j=0
do i=1,N_det do i=1,N_det
if (pruned(i)) cycle ! Pruned determinants if (pruned(i)) cycle ! Pruned determinants
j+=1 j+=1
ASSERT (sum(popcnt(psi_det(:,1,i))) == elec_alpha_num) ASSERT (sum(popcnt(psi_det(:,1,i))) == elec_alpha_num)
ASSERT (sum(popcnt(psi_det(:,2,i))) == elec_beta_num) ASSERT (sum(popcnt(psi_det(:,2,i))) == elec_beta_num)
buffer_det(:,:,j) = psi_det(:,:,i) buffer_det(:,:,j) = psi_det(:,:,i)
enddo enddo
N_det_old = j N_det_old = j
! Backup coefficients ! Backup coefficients
do k=1,N_states do k=1,N_states
j=0 j=0
do i=1,N_det do i=1,N_det
if (pruned(i)) cycle ! Pruned determinants if (pruned(i)) cycle ! Pruned determinants
j += 1 j += 1
buffer_coef(j,k) = psi_coef(i,k) buffer_coef(j,k) = psi_coef(i,k)
enddo enddo
ASSERT ( j == N_det_old ) ASSERT ( j == N_det_old )
enddo enddo
endif
! Update N_det ! Update N_det
N_det = N_det_old N_det = N_det_old
@ -164,17 +173,19 @@ subroutine copy_H_apply_buffer_to_wf
TOUCH psi_det_size TOUCH psi_det_size
endif endif
! Restore backup in resized array if (N_det_old > 0) then
do i=1,N_det_old ! Restore backup in resized array
psi_det(:,:,i) = buffer_det(:,:,i)
ASSERT (sum(popcnt(psi_det(:,1,i))) == elec_alpha_num)
ASSERT (sum(popcnt(psi_det(:,2,i))) == elec_beta_num )
enddo
do k=1,N_states
do i=1,N_det_old do i=1,N_det_old
psi_coef(i,k) = buffer_coef(i,k) psi_det(:,:,i) = buffer_det(:,:,i)
ASSERT (sum(popcnt(psi_det(:,1,i))) == elec_alpha_num)
ASSERT (sum(popcnt(psi_det(:,2,i))) == elec_beta_num )
enddo enddo
enddo do k=1,N_states
do i=1,N_det_old
psi_coef(i,k) = buffer_coef(i,k)
enddo
enddo
endif
! Copy new buffers ! Copy new buffers
@ -339,3 +350,33 @@ subroutine fill_H_apply_buffer_no_selection(n_selected,det_buffer,Nint,iproc)
call omp_unset_lock(H_apply_buffer_lock(1,iproc)) call omp_unset_lock(H_apply_buffer_lock(1,iproc))
end end
subroutine replace_wf(N_det_new, LDA, psi_coef_new, psi_det_new)
use omp_lib
implicit none
BEGIN_DOC
! Replaces the wave function.
! After calling this subroutine, N_det, psi_det and psi_coef need to be touched
END_DOC
integer, intent(in) :: N_det_new, LDA
double precision, intent(in) :: psi_coef_new(LDA,N_states)
integer(bit_kind), intent(in) :: psi_det_new(N_int,2,N_det_new)
integer :: i,j
PROVIDE H_apply_buffer_allocated
if (N_det_new <= 0) call abort() !irp_here//': N_det_new <= 0')
if (N_int <= 0) call abort() !irp_here//': N_int <= 0')
if (LDA < N_det_new) call abort() !irp_here//': LDA < N_det_new')
do j=0,nproc-1
H_apply_buffer(j)%N_det = 0
enddo
N_det = 0
SOFT_TOUCH N_det
call fill_H_apply_buffer_no_selection(N_det_new,psi_det_new,N_int,0)
call copy_h_apply_buffer_to_wf
psi_coef(1:N_det_new,1:N_states) = psi_coef_new(1:N_det_new,1:N_states)
end

30
src/kohn_sham/print_mos.irp.f Normal file
View File

@ -0,0 +1,30 @@
program print_mos
implicit none
integer :: i,nx
double precision :: r(3), xmax, dx, accu
double precision, allocatable :: mos_array(:)
double precision:: alpha,envelop
allocate(mos_array(mo_num))
xmax = 5.d0
nx = 1000
dx=xmax/dble(nx)
r = 0.d0
alpha = 0.5d0
do i = 1, nx
call give_all_mos_at_r(r,mos_array)
accu = mos_array(3)**2+mos_array(4)**2+mos_array(5)**2
accu = dsqrt(accu)
envelop = (1.d0 - dexp(-alpha * r(3)**2))
write(33,'(100(F16.10,X))')r(3), mos_array(1), mos_array(2), accu, envelop
r(3) += dx
enddo
end
double precision function f_mu(x)
implicit none
double precision, intent(in) :: x
end

18
src/nuclei/write_pt_charges.py
View File

@ -21,7 +21,7 @@ def mv_in_ezfio(ezfio,tmp):
os.system(cmdmv) os.system(cmdmv)
# Getting the EZFIO ##Getting the EZFIO
EZFIO=sys.argv[1] EZFIO=sys.argv[1]
EZFIO=EZFIO.replace("/", "") EZFIO=EZFIO.replace("/", "")
print(EZFIO) print(EZFIO)
@ -66,8 +66,20 @@ zip_in_ezfio(EZFIO,tmp)
tmp="pts_charge_coord" tmp="pts_charge_coord"
fcoord = open(tmp,'w') fcoord = open(tmp,'w')
fcoord.write(" 2\n") fcoord.write(" 2\n")
fcoord.write(" "+str(n_charges)+' 3\n') if(n_charges < 10):
#fcoord.write(" "+' 3 '+str(n_charges)+' \n') fcoord.write(" "+str(n_charges)+' 3\n')
elif(n_charges <100):
fcoord.write(" "+str(n_charges)+' 3\n')
elif(n_charges <1000):
fcoord.write(" "+str(n_charges)+' 3\n')
elif(n_charges <10000):
fcoord.write(" "+str(n_charges)+' 3\n')
elif(n_charges <100000):
fcoord.write(" "+str(n_charges)+' 3\n')
elif(n_charges <1000000):
fcoord.write(" "+str(n_charges)+' 3\n')
elif(n_charges <10000000):
fcoord.write(" "+str(n_charges)+' 3\n')
for i in range(n_charges): for i in range(n_charges):
fcoord.write(' '+coord_x[i]+'\n') fcoord.write(' '+coord_x[i]+'\n')
for i in range(n_charges): for i in range(n_charges):

549
src/tc_bi_ortho/dav_h_tc_s2.irp.f Normal file
View File

@ -0,0 +1,549 @@
! ---
subroutine davidson_hs2_nonsym_b1space(u_in, H_jj, s2_out,energies, sze, N_st, N_st_diag_in, converged, hcalc)
use mmap_module
BEGIN_DOC
! Generic modified-Davidson diagonalization
!
! H_jj : specific diagonal H matrix elements to diagonalize de Davidson
!
! u_in : guess coefficients on the various states. Overwritten on exit by right eigenvectors
!
! sze : Number of determinants
!
! N_st : Number of eigenstates
!
! N_st_diag_in : Number of states in which H is diagonalized. Assumed > N_st
!
! Initial guess vectors are not necessarily orthonormal
!
! hcalc subroutine to compute W = H U (see routine hcalc_template for template of input/output)
END_DOC
implicit none
integer, intent(in) :: sze, N_st, N_st_diag_in
double precision, intent(in) :: H_jj(sze)
logical, intent(inout) :: converged
double precision, intent(inout) :: u_in(sze,N_st_diag_in)
double precision, intent(out) :: energies(N_st)
double precision, intent(inout) :: s2_out(N_st)
external hcalc
character*(16384) :: write_buffer
integer :: iter, N_st_diag
integer :: i, j, k, l, m
integer :: iter2, itertot
logical :: disk_based
integer :: shift, shift2, itermax
integer :: nproc_target
integer :: order(N_st_diag_in)
double precision :: to_print(3,N_st)
double precision :: r1, r2, alpha
double precision :: cpu, wall
double precision :: cmax
double precision :: energy_shift(N_st_diag_in*davidson_sze_max)
double precision, allocatable :: U(:,:)
double precision, allocatable :: y(:,:), h(:,:), lambda(:), h_p(:,:), s2(:)
real, allocatable :: y_s(:,:)
double precision, allocatable :: s_(:,:), s_tmp(:,:)
double precision, allocatable :: residual_norm(:)
double precision :: lambda_tmp
integer, allocatable :: i_omax(:)
double precision, allocatable :: U_tmp(:), overlap(:), S_d(:,:)
double precision, allocatable :: W(:,:)
real, pointer :: S(:,:)
!double precision, pointer :: W(:,:)
double precision, external :: u_dot_v, u_dot_u
include 'constants.include.F'
N_st_diag = N_st_diag_in
! print*,'trial vector'
do i = 1, sze
if(isnan(u_in(i,1)))then
print*,'pb in input vector of davidson_general_ext_rout_nonsym_b1space'
print*,i,u_in(i,1)
stop
else if (dabs(u_in(i,1)).lt.1.d-16)then
u_in(i,1) = 0.d0
endif
enddo
!DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: U, W, S, y, y_s, S_d, h, lambda
if(N_st_diag*3 > sze) then
print *, 'error in Davidson :'
print *, 'Increase n_det_max_full to ', N_st_diag*3
stop -1
endif
itermax = max(2, min(davidson_sze_max, sze/N_st_diag)) + 1
provide threshold_nonsym_davidson
call write_time(6)
write(6,'(A)') ''
write(6,'(A)') 'Davidson Diagonalization'
write(6,'(A)') '------------------------'
write(6,'(A)') ''
! Find max number of cores to fit in memory
! -----------------------------------------
nproc_target = nproc
double precision :: rss
integer :: maxab
maxab = sze
m=1
disk_based = .False.
call resident_memory(rss)
do
r1 = 8.d0 * &! bytes
( dble(sze)*(N_st_diag*itermax) &! U
+ 1.5d0*dble(sze*m)*(N_st_diag*itermax) &! W, S
+ 4.5d0*(N_st_diag*itermax)**2 &! h,y,y_s,s_, s_tmp
+ 2.d0*(N_st_diag*itermax) &! s2,lambda
+ 1.d0*(N_st_diag) &! residual_norm
! In H_S2_u_0_nstates_zmq
+ 3.d0*(N_st_diag*N_det) &! u_t, v_t, s_t on collector
+ 3.d0*(N_st_diag*N_det) &! u_t, v_t, s_t on slave
+ 0.5d0*maxab &! idx0 in H_S2_u_0_nstates_openmp_work_*
+ nproc_target * &! In OMP section
( 1.d0*(N_int*maxab) &! buffer
+ 3.5d0*(maxab) ) &! singles_a, singles_b, doubles, idx
) / 1024.d0**3
if(nproc_target == 0) then
call check_mem(r1, irp_here)
nproc_target = 1
exit
endif
if(r1+rss < qp_max_mem) then
exit
endif
if(itermax > 4) then
itermax = itermax - 1
! else if (m==1.and.disk_based_davidson) then
! m = 0
! disk_based = .True.
! itermax = 6
else
nproc_target = nproc_target - 1
endif
enddo
nthreads_davidson = nproc_target
TOUCH nthreads_davidson
call write_int(6, N_st, 'Number of states')
call write_int(6, N_st_diag, 'Number of states in diagonalization')
call write_int(6, sze, 'Number of basis functions')
call write_int(6, nproc_target, 'Number of threads for diagonalization')
call write_double(6, r1, 'Memory(Gb)')
if(disk_based) then
print *, 'Using swap space to reduce RAM'
endif
!---------------
write(6,'(A)') ''
write_buffer = '====='
do i=1,N_st
write_buffer = trim(write_buffer)//' ================ =========== ==========='
enddo
write(6,'(A)') write_buffer(1:6+41*N_st)
write_buffer = 'Iter'
do i=1,N_st
write_buffer = trim(write_buffer)//' Energy S^2 Residual '
enddo
write(6,'(A)') write_buffer(1:6+41*N_st)
write_buffer = '====='
do i=1,N_st
write_buffer = trim(write_buffer)//' ================ =========== ==========='
enddo
write(6,'(A)') write_buffer(1:6+41*N_st)
! ---
allocate( W(sze,N_st_diag*itermax), S(sze,N_st_diag*itermax) )
allocate( &
! Large
U(sze,N_st_diag*itermax), &
S_d(sze,N_st_diag), &
! Small
h(N_st_diag*itermax,N_st_diag*itermax), &
h_p(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_tmp(N_st_diag*itermax,N_st_diag*itermax), &
lambda(N_st_diag*itermax), &
residual_norm(N_st_diag), &
i_omax(N_st), &
s2(N_st_diag*itermax), &
y_s(N_st_diag*itermax,N_st_diag*itermax) &
)
U = 0.d0
h = 0.d0
y = 0.d0
s_ = 0.d0
s_tmp = 0.d0
lambda = 0.d0
residual_norm = 0.d0
ASSERT (N_st > 0)
ASSERT (N_st_diag >= N_st)
ASSERT (sze > 0)
! Davidson iterations
! ===================
converged = .False.
! Initialize from N_st to N_st_diag with gaussian random numbers
! to be sure to have overlap with any eigenvectors
do k = N_st+1, N_st_diag
u_in(k,k) = 10.d0
do i = 1, sze
call random_number(r1)
call random_number(r2)
r1 = dsqrt(-2.d0*dlog(r1))
r2 = dtwo_pi*r2
u_in(i,k) = r1*dcos(r2)
enddo
enddo
! Normalize all states
do k = 1, N_st_diag
call normalize(u_in(1,k), sze)
enddo
! Copy from the guess input "u_in" to the working vectors "U"
do k = 1, N_st_diag
do i = 1, sze
U(i,k) = u_in(i,k)
enddo
enddo
! ---
itertot = 0
do while (.not.converged)
itertot = itertot + 1
if(itertot == 8) then
exit
endif
do iter = 1, itermax-1
shift = N_st_diag * (iter-1)
shift2 = N_st_diag * iter
if( (iter > 1) .or. (itertot == 1) ) then
! Gram-Schmidt to orthogonalize all new guess with the previous vectors
call ortho_qr(U, size(U, 1), sze, shift2)
call ortho_qr(U, size(U, 1), sze, shift2)
! W = H U
! call hcalc(W(1,shift+1), U(1,shift+1), N_st_diag, sze)
call hcalc(W(1,shift+1),S_d,U(1,shift+1),N_st_diag,sze)
S(1:sze,shift+1:shift+N_st_diag) = real(S_d(1:sze,1:N_st_diag))
else
! Already computed in update below
continue
endif
! Compute s_kl = <u_k | S_l> = <u_k| S2 |u_l>
! -------------------------------------------
!$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(i,j,k) COLLAPSE(2)
do j=1,shift2
do i=1,shift2
s_(i,j) = 0.d0
do k=1,sze
s_(i,j) = s_(i,j) + U(k,i) * dble(S(k,j))
enddo
enddo
enddo
!$OMP END PARALLEL DO
! Compute h_kl = <u_k | W_l> = <u_k| H |u_l>
! -------------------------------------------
call dgemm( 'T', 'N', shift2, shift2, sze, 1.d0 &
, U, size(U, 1), W, size(W, 1) &
, 0.d0, h, size(h, 1) )
! Penalty method
! --------------
if (s2_eig) then
h_p = s_
do k=1,shift2
h_p(k,k) = h_p(k,k) - expected_s2
enddo
if (only_expected_s2) then
alpha = 0.1d0
h_p = h + alpha*h_p
else
alpha = 0.0001d0
h_p = h + alpha*h_p
endif
else
h_p = h
alpha = 0.d0
endif
! Diagonalize h y = lambda y
! ---------------------------
call diag_nonsym_right(shift2, h_p(1,1), size(h_p, 1), y(1,1), size(y, 1), lambda(1), size(lambda, 1))
do k = 1, N_st_diag
! print*,'lambda(k) before = ',lambda(k)
lambda(k) = 0.d0
do l = 1, shift2
do m = 1, shift2
lambda(k) += y(m,k) * h(m,l) * y(l,k)
enddo
enddo
! print*,'lambda(k) new = ',lambda(k)
enddo
! Compute S2 for each eigenvector
! -------------------------------
call dgemm('N','N',shift2,shift2,shift2, &
1.d0, s_, size(s_,1), y, size(y,1), &
0.d0, s_tmp, size(s_tmp,1))
call dgemm('T','N',shift2,shift2,shift2, &
1.d0, y, size(y,1), s_tmp, size(s_tmp,1), &
0.d0, s_, size(s_,1))
do k=1,shift2
s2(k) = s_(k,k)
enddo
! Express eigenvectors of h in the determinant basis:
! ---------------------------------------------------
! y(:,k) = rk
! U(:,k) = Bk
! U(:,shift2+k) = Rk = Bk x rk
call dgemm( 'N', 'N', sze, N_st_diag, shift2, 1.d0 &
, U, size(U, 1), y, size(y, 1) &
, 0.d0, U(1,shift2+1), size(U, 1) )
do k = 1, N_st_diag
call normalize(U(1,shift2+k), sze)
enddo
! ---
! select the max overlap
!
! start test ------------------------------------------------------------------------
!
!double precision, allocatable :: Utest(:,:), Otest(:)
!allocate( Utest(sze,shift2), Otest(shift2) )
!call dgemm( 'N', 'N', sze, shift2, shift2, 1.d0 &
! , U, size(U, 1), y, size(y, 1), 0.d0, Utest(1,1), size(Utest, 1) )
!do k = 1, shift2
! call normalize(Utest(1,k), sze)
!enddo
!do j = 1, sze
! write(455, '(100(1X, F16.10))') (Utest(j,k), k=1,shift2)
!enddo
!do k = 1, shift2
! Otest(k) = 0.d0
! do i = 1, sze
! Otest(k) += Utest(i,k) * u_in(i,1)
! enddo
! Otest(k) = dabs(Otest(k))
! print *, ' Otest =', k, Otest(k), lambda(k)
!enddo
!deallocate(Utest, Otest)
!
! end test ------------------------------------------------------------------------
!
! TODO
! state_following is more efficient
do l = 1, N_st
allocate( overlap(N_st_diag) )
do k = 1, N_st_diag
overlap(k) = 0.d0
do i = 1, sze
overlap(k) = overlap(k) + U(i,shift2+k) * u_in(i,l)
enddo
overlap(k) = dabs(overlap(k))
!print *, ' overlap =', k, overlap(k)
enddo
lambda_tmp = 0.d0
do k = 1, N_st_diag
if(overlap(k) .gt. lambda_tmp) then
i_omax(l) = k
lambda_tmp = overlap(k)
endif
enddo
deallocate(overlap)
if(lambda_tmp .lt. 0.7d0) then
print *, ' very small overlap ...', l, i_omax(l)
print *, ' max overlap = ', lambda_tmp
stop
endif
if(i_omax(l) .ne. l) then
print *, ' !!! WARNONG !!!'
print *, ' index of state', l, i_omax(l)
endif
enddo
! y(:,k) = rk
! W(:,k) = H x Bk
! W(:,shift2+k) = H x Bk x rk
! = Wk
call dgemm( 'N', 'N', sze, N_st_diag, shift2, 1.d0 &
, W, size(W, 1), y, size(y, 1) &
, 0.d0, W(1,shift2+1), size(W, 1) )
! ---
! Compute residual vector and davidson step
! -----------------------------------------
!$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(i,k)
do k = 1, N_st_diag
do i = 1, sze
U(i,shift2+k) = (lambda(k) * U(i,shift2+k) - W(i,shift2+k)) / max(H_jj(i)-lambda(k), 1.d-2)
enddo
if(k <= N_st) then
l = k
residual_norm(k) = u_dot_u(U(1,shift2+l), sze)
to_print(1,k) = lambda(l)
to_print(2,k) = s2(l)
to_print(3,k) = residual_norm(l)
endif
enddo
!$OMP END PARALLEL DO
!residual_norm(1) = u_dot_u(U(1,shift2+1), sze)
!to_print(1,1) = lambda(1)
!to_print(2,1) = residual_norm(1)
if( (itertot > 1) .and. (iter == 1) ) then
!don't print
continue
else
write(*, '(1X, I3, 1X, 100(1X, F16.10, 1X, F16.10, 1X, F16.10))') iter-1, to_print(1:3,1:N_st)
endif
! Check convergence
if(iter > 1) then
converged = dabs(maxval(residual_norm(1:N_st))) < threshold_nonsym_davidson
endif
do k = 1, N_st
if(residual_norm(k) > 1.e8) then
print *, 'Davidson failed'
stop -1
endif
enddo
if(converged) then
exit
endif
logical, external :: qp_stop
if(qp_stop()) then
converged = .True.
exit
endif
enddo ! loop over iter
! Re-contract U and update W
! --------------------------------
call dgemm( 'N', 'N', sze, N_st_diag, shift2, 1.d0 &
, W, size(W, 1), y, size(y, 1) &
, 0.d0, u_in, size(u_in, 1) )
do k = 1, N_st_diag
do i = 1, sze
W(i,k) = u_in(i,k)
enddo
enddo
call dgemm( 'N', 'N', sze, N_st_diag, shift2, 1.d0 &
, U, size(U, 1), y, size(y, 1) &
, 0.d0, u_in, size(u_in, 1) )
do k = 1, N_st_diag
do i = 1, sze
U(i,k) = u_in(i,k)
enddo
enddo
call ortho_qr(U, size(U, 1), sze, N_st_diag)
call ortho_qr(U, size(U, 1), sze, N_st_diag)
do j = 1, N_st_diag
k = 1
do while( (k < sze) .and. (U(k,j) == 0.d0) )
k = k+1
enddo
if(U(k,j) * u_in(k,j) < 0.d0) then
do i = 1, sze
W(i,j) = -W(i,j)
enddo
endif
enddo
enddo ! loop over while
! ---
do k = 1, N_st
energies(k) = lambda(k)
s2_out(k) = s2(k)
enddo
write_buffer = '====='
do i = 1, N_st
write_buffer = trim(write_buffer)//' ================ ==========='
enddo
write(6,'(A)') trim(write_buffer)
write(6,'(A)') ''
call write_time(6)
deallocate(W)
deallocate(U, h, y, lambda, residual_norm, i_omax)
FREE nthreads_davidson
end subroutine davidson_general_ext_rout_nonsym_b1space
! ---

769
src/tc_bi_ortho/h_tc_s2_u0.irp.f Normal file
View File

@ -0,0 +1,769 @@
subroutine get_H_tc_s2_l0_r0(l_0,r_0,N_st,sze,energies, s2)
use bitmasks
implicit none
BEGIN_DOC
! Computes $e_0 = \langle l_0 | H | r_0\rangle$.
!
! Computes $s_0 = \langle l_0 | S^2 | r_0\rangle$.
!
! Assumes that the determinants are in psi_det
!
! istart, iend, ishift, istep are used in ZMQ parallelization.
END_DOC
integer, intent(in) :: N_st,sze
double precision, intent(in) :: l_0(sze,N_st), r_0(sze,N_st)
double precision, intent(out) :: energies(N_st), s2(N_st)
logical :: do_right
integer :: istate
double precision, allocatable :: s_0(:,:), v_0(:,:)
double precision :: u_dot_v, norm
allocate(s_0(sze,N_st), v_0(sze,N_st))
do_right = .True.
call H_tc_s2_u_0_opt(v_0,s_0,r_0,N_st,sze)
do istate = 1, N_st
norm = u_dot_v(l_0(1,istate),r_0(1,istate),sze)
energies(istate) = u_dot_v(l_0(1,istate),v_0(1,istate),sze)/norm
s2(istate) = u_dot_v(l_0(1,istate),s_0(1,istate),sze)/norm
enddo
end
subroutine H_tc_s2_u_0_opt(v_0,s_0,u_0,N_st,sze)
use bitmasks
implicit none
BEGIN_DOC
! Computes $v_0 = H | u_0\rangle$.
!
! Assumes that the determinants are in psi_det
!
! istart, iend, ishift, istep are used in ZMQ parallelization.
END_DOC
integer, intent(in) :: N_st,sze
double precision, intent(inout) :: v_0(sze,N_st), u_0(sze,N_st), s_0(sze,N_st)
logical :: do_right
do_right = .True.
call H_tc_s2_u_0_nstates_openmp(v_0,s_0,u_0,N_st,sze, do_right)
end
subroutine H_tc_s2_dagger_u_0_opt(v_0,s_0,u_0,N_st,sze)
use bitmasks
implicit none
BEGIN_DOC
! Computes $v_0 = H | u_0\rangle$.
!
! Assumes that the determinants are in psi_det
!
! istart, iend, ishift, istep are used in ZMQ parallelization.
END_DOC
integer, intent(in) :: N_st,sze
double precision, intent(inout) :: v_0(sze,N_st), u_0(sze,N_st), s_0(sze,N_st)
logical :: do_right
do_right = .False.
call H_tc_s2_u_0_nstates_openmp(v_0,s_0,u_0,N_st,sze, do_right)
end
subroutine H_tc_s2_u_0_nstates_openmp(v_0,s_0,u_0,N_st,sze, do_right)
use bitmasks
implicit none
BEGIN_DOC
! Computes $v_0 = H | u_0\rangle$.
!
! Assumes that the determinants are in psi_det
!
! istart, iend, ishift, istep are used in ZMQ parallelization.
!
! if do_right == True then you compute H_TC |Psi>, else H_TC^T |Psi>
END_DOC
integer, intent(in) :: N_st,sze
double precision, intent(inout) :: v_0(sze,N_st), u_0(sze,N_st), s_0(sze,N_st)
logical, intent(in) :: do_right
integer :: k
double precision, allocatable :: u_t(:,:), v_t(:,:), s_t(:,:)
!DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: u_t
allocate(u_t(N_st,N_det),v_t(N_st,N_det),s_t(N_st,N_det))
do k=1,N_st
call dset_order(u_0(1,k),psi_bilinear_matrix_order,N_det)
enddo
v_t = 0.d0
s_t = 0.d0
call dtranspose( &
u_0, &
size(u_0, 1), &
u_t, &
size(u_t, 1), &
N_det, N_st)
call H_tc_s2_u_0_nstates_openmp_work(v_t,s_t,u_t,N_st,sze,1,N_det,0,1, do_right)
deallocate(u_t)
call dtranspose( &
v_t, &
size(v_t, 1), &
v_0, &
size(v_0, 1), &
N_st, N_det)
call dtranspose( &
s_t, &
size(s_t, 1), &
s_0, &
size(s_0, 1), &
N_st, N_det)
deallocate(v_t,s_t)
do k=1,N_st
call dset_order(v_0(1,k),psi_bilinear_matrix_order_reverse,N_det)
call dset_order(s_0(1,k),psi_bilinear_matrix_order_reverse,N_det)
call dset_order(u_0(1,k),psi_bilinear_matrix_order_reverse,N_det)
enddo
end
subroutine H_tc_s2_u_0_nstates_openmp_work(v_t,s_t,u_t,N_st,sze,istart,iend,ishift,istep, do_right)
use bitmasks
implicit none
BEGIN_DOC
! Computes $v_t = H | u_t\rangle$
!
! Default should be 1,N_det,0,1
!
! if do_right == True then you compute H_TC |Psi>, else H_TC^T |Psi>
END_DOC
integer, intent(in) :: N_st,sze,istart,iend,ishift,istep
double precision, intent(in) :: u_t(N_st,N_det)
logical, intent(in) :: do_right
double precision, intent(out) :: v_t(N_st,sze), s_t(N_st,sze)
PROVIDE ref_bitmask_energy N_int
select case (N_int)
case (1)
call H_tc_s2_u_0_nstates_openmp_work_1(v_t,s_t,u_t,N_st,sze,istart,iend,ishift,istep,do_right)
case (2)
call H_tc_s2_u_0_nstates_openmp_work_2(v_t,s_t,u_t,N_st,sze,istart,iend,ishift,istep,do_right)
case (3)
call H_tc_s2_u_0_nstates_openmp_work_3(v_t,s_t,u_t,N_st,sze,istart,iend,ishift,istep,do_right)
case (4)
call H_tc_s2_u_0_nstates_openmp_work_4(v_t,s_t,u_t,N_st,sze,istart,iend,ishift,istep,do_right)
case default
call H_tc_s2_u_0_nstates_openmp_work_N_int(v_t,s_t,u_t,N_st,sze,istart,iend,ishift,istep,do_right)
end select
end
BEGIN_TEMPLATE
subroutine H_tc_s2_u_0_nstates_openmp_work_$N_int(v_t,s_t,u_t,N_st,sze,istart,iend,ishift,istep,do_right)
use bitmasks
implicit none
BEGIN_DOC
! Computes $v_t = H | u_t \\rangle$ and $s_t = S^2 | u_t\\rangle$
!
! Default should be 1,N_det,0,1
!
! if do_right == True then you compute H_TC |Psi>, else H_TC^T |Psi>
END_DOC
integer, intent(in) :: N_st,sze,istart,iend,ishift,istep
double precision, intent(in) :: u_t(N_st,N_det)
logical, intent(in) :: do_right
double precision, intent(out) :: v_t(N_st,sze), s_t(N_st,sze)
double precision :: hij, sij
integer :: i,j,k,l,kk
integer :: k_a, k_b, l_a, l_b, m_a, m_b
integer :: istate
integer :: krow, kcol, krow_b, kcol_b
integer :: lrow, lcol
integer :: mrow, mcol
integer(bit_kind) :: spindet($N_int)
integer(bit_kind) :: tmp_det($N_int,2)
integer(bit_kind) :: tmp_det2($N_int,2)
integer(bit_kind) :: tmp_det3($N_int,2)
integer(bit_kind), allocatable :: buffer(:,:)
integer :: n_doubles
integer, allocatable :: doubles(:)
integer, allocatable :: singles_a(:)
integer, allocatable :: singles_b(:)
integer, allocatable :: idx(:), idx0(:)
integer :: maxab, n_singles_a, n_singles_b, kcol_prev
integer*8 :: k8
logical :: compute_singles
integer*8 :: last_found, left, right, right_max
double precision :: rss, mem, ratio
double precision, allocatable :: utl(:,:)
integer, parameter :: block_size=128
logical :: u_is_sparse
! call resident_memory(rss)
! mem = dble(singles_beta_csc_size) / 1024.d0**3
!
! compute_singles = (mem+rss > qp_max_mem)
!
! if (.not.compute_singles) then
! provide singles_beta_csc
! endif
compute_singles=.True.
maxab = max(N_det_alpha_unique, N_det_beta_unique)+1
allocate(idx0(maxab))
do i=1,maxab
idx0(i) = i
enddo
! Prepare the array of all alpha single excitations
! -------------------------------------------------
PROVIDE N_int nthreads_davidson
!$OMP PARALLEL DEFAULT(SHARED) NUM_THREADS(nthreads_davidson) &
!$OMP SHARED(psi_bilinear_matrix_rows, N_det, &
!$OMP psi_bilinear_matrix_columns, &
!$OMP psi_det_alpha_unique, psi_det_beta_unique, &
!$OMP n_det_alpha_unique, n_det_beta_unique, N_int, &
!$OMP psi_bilinear_matrix_transp_rows, &
!$OMP psi_bilinear_matrix_transp_columns, &
!$OMP psi_bilinear_matrix_transp_order, N_st, &
!$OMP psi_bilinear_matrix_order_transp_reverse, &
!$OMP psi_bilinear_matrix_columns_loc, &
!$OMP psi_bilinear_matrix_transp_rows_loc, &
!$OMP istart, iend, istep, irp_here, v_t, s_t, &
!$OMP ishift, idx0, u_t, maxab, compute_singles, &
!$OMP singles_alpha_csc,singles_alpha_csc_idx, &
!$OMP singles_beta_csc,singles_beta_csc_idx) &
!$OMP PRIVATE(krow, kcol, tmp_det, spindet, k_a, k_b, i, &
!$OMP lcol, lrow, l_a, l_b, utl, kk, u_is_sparse, &
!$OMP buffer, doubles, n_doubles, umax, &
!$OMP tmp_det2, hij, sij, idx, l, kcol_prev,hmono, htwoe, hthree, &
!$OMP singles_a, n_singles_a, singles_b, ratio, &
!$OMP n_singles_b, k8, last_found,left,right,right_max)
! Alpha/Beta double excitations
! =============================
allocate( buffer($N_int,maxab), &
singles_a(maxab), &
singles_b(maxab), &
doubles(maxab), &
idx(maxab), utl(N_st,block_size))
kcol_prev=-1
! Check if u has multiple zeros
kk=1 ! Avoid division by zero
!$OMP DO
do k=1,N_det
umax = 0.d0
do l=1,N_st
umax = max(umax, dabs(u_t(l,k)))
enddo
if (umax < 1.d-20) then
!$OMP ATOMIC
kk = kk+1
endif
enddo
!$OMP END DO
u_is_sparse = N_det / kk < 20 ! 5%
ASSERT (iend <= N_det)
ASSERT (istart > 0)
ASSERT (istep > 0)
!$OMP DO SCHEDULE(guided,64)
do k_a=istart+ishift,iend,istep ! Loop over all determinants (/!\ not in psidet order)
krow = psi_bilinear_matrix_rows(k_a) ! Index of alpha part of determinant k_a
ASSERT (krow <= N_det_alpha_unique)
kcol = psi_bilinear_matrix_columns(k_a) ! Index of beta part of determinant k_a
ASSERT (kcol <= N_det_beta_unique)
tmp_det(1:$N_int,1) = psi_det_alpha_unique(1:$N_int, krow)
if (kcol /= kcol_prev) then
tmp_det(1:$N_int,2) = psi_det_beta_unique (1:$N_int, kcol)
if (compute_singles) then
call get_all_spin_singles_$N_int( &
psi_det_beta_unique, idx0, &
tmp_det(1,2), N_det_beta_unique, &
singles_b, n_singles_b)
else
n_singles_b = 0
!DIR$ LOOP COUNT avg(1000)
do k8=singles_beta_csc_idx(kcol),singles_beta_csc_idx(kcol+1)-1
n_singles_b = n_singles_b+1
singles_b(n_singles_b) = singles_beta_csc(k8)
enddo
endif
endif
kcol_prev = kcol
! -> Here, tmp_det is determinant k_a
! Loop over singly excited beta columns
! -------------------------------------
!DIR$ LOOP COUNT avg(1000)
do i=1,n_singles_b
lcol = singles_b(i)
tmp_det2(1:$N_int,2) = psi_det_beta_unique(1:$N_int, lcol)
! tmp_det2 is a single excitation of tmp_det in the beta spin
! the alpha part is not defined yet
!---
! if (compute_singles) then
l_a = psi_bilinear_matrix_columns_loc(lcol)
ASSERT (l_a <= N_det)
! rows : | 1 2 3 4 | 1 3 4 6 | .... | 1 2 4 5 |
! cols : | 1 1 1 1 | 2 2 2 2 | .... | 8 8 8 8 |
! index : | 1 2 3 4 | 5 6 7 8 | .... | 58 59 60 61 |
! ^ ^
! | |
! l_a N_det
! l_a is the index in the big vector os size Ndet of the position of the first element of column lcol
! Below we identify all the determinants with the same beta part
!DIR$ UNROLL(8)
!DIR$ LOOP COUNT avg(50000)
do j=1,psi_bilinear_matrix_columns_loc(lcol+1) - psi_bilinear_matrix_columns_loc(lcol)
lrow = psi_bilinear_matrix_rows(l_a)
ASSERT (lrow <= N_det_alpha_unique)
buffer(1:$N_int,j) = psi_det_alpha_unique(1:$N_int, lrow) ! hot spot
ASSERT (l_a <= N_det)
idx(j) = l_a
l_a = l_a+1
enddo
j = j-1
! Get all single excitations from tmp_det(1,1) to buffer(1,?)
call get_all_spin_singles_$N_int( &
buffer, idx, tmp_det(1,1), j, &
singles_a, n_singles_a )
! Loop over alpha singles
! -----------------------
double precision :: umax
!DIR$ LOOP COUNT avg(1000)
do k = 1,n_singles_a,block_size
umax = 0.d0
! Prefetch u_t(:,l_a)
if (u_is_sparse) then
do kk=0,block_size-1
if (k+kk > n_singles_a) exit
l_a = singles_a(k+kk)
ASSERT (l_a <= N_det)
do l=1,N_st
utl(l,kk+1) = u_t(l,l_a)
umax = max(umax, dabs(utl(l,kk+1)))
enddo
enddo
else
do kk=0,block_size-1
if (k+kk > n_singles_a) exit
l_a = singles_a(k+kk)
ASSERT (l_a <= N_det)
utl(:,kk+1) = u_t(:,l_a)
enddo
umax = 1.d0
endif
if (umax < 1.d-20) cycle
do kk=0,block_size-1
if (k+kk > n_singles_a) exit
l_a = singles_a(k+kk)
lrow = psi_bilinear_matrix_rows(l_a)
ASSERT (lrow <= N_det_alpha_unique)
tmp_det2(1:$N_int,1) = psi_det_alpha_unique(1:$N_int, lrow)
! call i_H_j( tmp_det, tmp_det2, $N_int, hij) ! double alpha-beta
if(do_right)then
call htilde_mu_mat_opt_bi_ortho_tot(tmp_det,tmp_det2,$N_int,hij)
else
call htilde_mu_mat_opt_bi_ortho_tot(tmp_det2,tmp_det,$N_int,hij)
endif
call get_s2(tmp_det,tmp_det2,$N_int,sij)
!DIR$ LOOP COUNT AVG(4)
do l=1,N_st
v_t(l,k_a) = v_t(l,k_a) + hij * utl(l,kk+1)
s_t(l,k_a) = s_t(l,k_a) + sij * utl(l,kk+1)
enddo
enddo
enddo
enddo
enddo
!$OMP END DO
!$OMP DO SCHEDULE(guided,64)
do k_a=istart+ishift,iend,istep
! Single and double alpha excitations
! ===================================
! Initial determinant is at k_a in alpha-major representation
! -----------------------------------------------------------------------
krow = psi_bilinear_matrix_rows(k_a)
ASSERT (krow <= N_det_alpha_unique)
kcol = psi_bilinear_matrix_columns(k_a)
ASSERT (kcol <= N_det_beta_unique)
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)
! Initial determinant is at k_b in beta-major representation
! ----------------------------------------------------------------------
k_b = psi_bilinear_matrix_order_transp_reverse(k_a)
ASSERT (k_b <= N_det)
spindet(1:$N_int) = tmp_det(1:$N_int,1)
! Loop inside the beta column to gather all the connected alphas
lcol = psi_bilinear_matrix_columns(k_a)
l_a = psi_bilinear_matrix_columns_loc(lcol)
!DIR$ LOOP COUNT avg(200000)
do i=1,N_det_alpha_unique
if (l_a > N_det) exit
lcol = psi_bilinear_matrix_columns(l_a)
if (lcol /= kcol) exit
lrow = psi_bilinear_matrix_rows(l_a)
ASSERT (lrow <= N_det_alpha_unique)
buffer(1:$N_int,i) = psi_det_alpha_unique(1:$N_int, lrow) ! Hot spot
idx(i) = l_a
l_a = l_a+1
enddo
i = i-1
call get_all_spin_singles_and_doubles_$N_int( &
buffer, idx, spindet, i, &
singles_a, doubles, n_singles_a, n_doubles )
! Compute Hij for all alpha singles
! ----------------------------------
tmp_det2(1:$N_int,2) = psi_det_beta_unique (1:$N_int, kcol)
!DIR$ LOOP COUNT avg(1000)
do i=1,n_singles_a,block_size
umax = 0.d0
! Prefetch u_t(:,l_a)
if (u_is_sparse) then
do kk=0,block_size-1
if (i+kk > n_singles_a) exit
l_a = singles_a(i+kk)
ASSERT (l_a <= N_det)
do l=1,N_st
utl(l,kk+1) = u_t(l,l_a)
umax = max(umax, dabs(utl(l,kk+1)))
enddo
enddo
else
do kk=0,block_size-1
if (i+kk > n_singles_a) exit
l_a = singles_a(i+kk)
ASSERT (l_a <= N_det)
utl(:,kk+1) = u_t(:,l_a)
enddo
umax = 1.d0
endif
if (umax < 1.d-20) cycle
do kk=0,block_size-1
if (i+kk > n_singles_a) exit
l_a = singles_a(i+kk)
lrow = psi_bilinear_matrix_rows(l_a)
ASSERT (lrow <= N_det_alpha_unique)
tmp_det2(1:$N_int,1) = psi_det_alpha_unique(1:$N_int, lrow)
! call i_h_j_single_spin( tmp_det, tmp_det2, $N_int, 1, hij)
if(do_right)then
call htilde_mu_mat_opt_bi_ortho_tot(tmp_det,tmp_det2,$N_int,hij)
else
call htilde_mu_mat_opt_bi_ortho_tot(tmp_det2,tmp_det,$N_int,hij)
endif
!DIR$ LOOP COUNT AVG(4)
do l=1,N_st
v_t(l,k_a) = v_t(l,k_a) + hij * utl(l,kk+1)
enddo
enddo
enddo
! Compute Hij for all alpha doubles
! ----------------------------------
!DIR$ LOOP COUNT avg(50000)
do i=1,n_doubles,block_size
umax = 0.d0
! Prefetch u_t(:,l_a)
if (u_is_sparse) then
do kk=0,block_size-1
if (i+kk > n_doubles) exit
l_a = doubles(i+kk)
ASSERT (l_a <= N_det)
do l=1,N_st
utl(l,kk+1) = u_t(l,l_a)
umax = max(umax, dabs(utl(l,kk+1)))
enddo
enddo
else
do kk=0,block_size-1
if (i+kk > n_doubles) exit
l_a = doubles(i+kk)
ASSERT (l_a <= N_det)
utl(:,kk+1) = u_t(:,l_a)
enddo
umax = 1.d0
endif
if (umax < 1.d-20) cycle
do kk=0,block_size-1
if (i+kk > n_doubles) exit
l_a = doubles(i+kk)
lrow = psi_bilinear_matrix_rows(l_a)
ASSERT (lrow <= N_det_alpha_unique)
tmp_det2(1:N_int,1) = psi_det_alpha_unique(1:N_int, lrow)
! call i_H_j( tmp_det, tmp_det2, $N_int, hij)
! call i_H_j_double_spin( tmp_det(1,1), psi_det_alpha_unique(1, lrow), $N_int, hij)
if(do_right)then
call htilde_mu_mat_opt_bi_ortho_tot(tmp_det,tmp_det2,$N_int,hij)
else
call htilde_mu_mat_opt_bi_ortho_tot(tmp_det2,tmp_det,$N_int,hij)
endif
!DIR$ LOOP COUNT AVG(4)
do l=1,N_st
v_t(l,k_a) = v_t(l,k_a) + hij * utl(l,kk+1)
enddo
enddo
enddo
! Single and double beta excitations
! ==================================
! Initial determinant is at k_a in alpha-major representation
! -----------------------------------------------------------------------
krow = psi_bilinear_matrix_rows(k_a)
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)
spindet(1:$N_int) = tmp_det(1:$N_int,2)
! Initial determinant is at k_b in beta-major representation
! -----------------------------------------------------------------------
k_b = psi_bilinear_matrix_order_transp_reverse(k_a)
ASSERT (k_b <= N_det)
! Loop inside the alpha row to gather all the connected betas
lrow = psi_bilinear_matrix_transp_rows(k_b)
l_b = psi_bilinear_matrix_transp_rows_loc(lrow)
!DIR$ LOOP COUNT avg(200000)
do i=1,N_det_beta_unique
if (l_b > N_det) exit
lrow = psi_bilinear_matrix_transp_rows(l_b)
if (lrow /= krow) exit
lcol = psi_bilinear_matrix_transp_columns(l_b)
ASSERT (lcol <= N_det_beta_unique)
buffer(1:$N_int,i) = psi_det_beta_unique(1:$N_int, lcol)
idx(i) = l_b
l_b = l_b+1
enddo
i = i-1
call get_all_spin_singles_and_doubles_$N_int( &
buffer, idx, spindet, i, &
singles_b, doubles, n_singles_b, n_doubles )
! Compute Hij for all beta singles
! ----------------------------------
tmp_det2(1:$N_int,1) = psi_det_alpha_unique(1:$N_int, krow)
!DIR$ LOOP COUNT avg(1000)
do i=1,n_singles_b,block_size
umax = 0.d0
if (u_is_sparse) then
do kk=0,block_size-1
if (i+kk > n_singles_b) exit
l_b = singles_b(i+kk)
l_a = psi_bilinear_matrix_transp_order(l_b)
ASSERT (l_b <= N_det)
ASSERT (l_a <= N_det)
do l=1,N_st
utl(l,kk+1) = u_t(l,l_a)
umax = max(umax, dabs(utl(l,kk+1)))
enddo
enddo
else
do kk=0,block_size-1
if (i+kk > n_singles_b) exit
l_b = singles_b(i+kk)
l_a = psi_bilinear_matrix_transp_order(l_b)
ASSERT (l_b <= N_det)
ASSERT (l_a <= N_det)
utl(:,kk+1) = u_t(:,l_a)
enddo
umax = 1.d0
endif
if (umax < 1.d-20) cycle
do kk=0,block_size-1
if (i+kk > n_singles_b) exit
l_b = singles_b(i+kk)
l_a = psi_bilinear_matrix_transp_order(l_b)
lcol = psi_bilinear_matrix_transp_columns(l_b)
ASSERT (lcol <= N_det_beta_unique)
tmp_det2(1:$N_int,2) = psi_det_beta_unique (1:$N_int, lcol)
! call i_H_j_single_spin( tmp_det, tmp_det2, $N_int, 2, hij)
if(do_right)then
call htilde_mu_mat_opt_bi_ortho_tot(tmp_det,tmp_det2,$N_int,hij)
else
call htilde_mu_mat_opt_bi_ortho_tot(tmp_det2,tmp_det,$N_int,hij)
endif
!DIR$ LOOP COUNT AVG(4)
do l=1,N_st
v_t(l,k_a) = v_t(l,k_a) + hij * utl(l,kk+1)
enddo
enddo
enddo
! Compute Hij for all beta doubles
! ----------------------------------
!DIR$ LOOP COUNT avg(50000)
do i=1,n_doubles,block_size
umax = 0.d0
if (u_is_sparse) then
do kk=0,block_size-1
if (i+kk > n_doubles) exit
l_b = doubles(i+kk)
l_a = psi_bilinear_matrix_transp_order(l_b)
ASSERT (l_b <= N_det)
ASSERT (l_a <= N_det)
do l=1,N_st
utl(l,kk+1) = u_t(l,l_a)
umax = max(umax, dabs(utl(l,kk+1)))
enddo
enddo
else
do kk=0,block_size-1
if (i+kk > n_doubles) exit
l_b = doubles(i+kk)
l_a = psi_bilinear_matrix_transp_order(l_b)
ASSERT (l_b <= N_det)
ASSERT (l_a <= N_det)
utl(:,kk+1) = u_t(:,l_a)
enddo
umax = 1.d0
endif
if (umax < 1.d-20) cycle
do kk=0,block_size-1
if (i+kk > n_doubles) exit
l_b = doubles(i+kk)
l_a = psi_bilinear_matrix_transp_order(l_b)
lcol = psi_bilinear_matrix_transp_columns(l_b)
ASSERT (lcol <= N_det_beta_unique)
tmp_det2(1:N_int,2) = psi_det_beta_unique(1:N_int, lcol)
! call i_H_j( tmp_det, tmp_det2, $N_int, hij)
! call i_H_j_double_spin( tmp_det(1,2), psi_det_beta_unique(1, lcol), $N_int, hij)
if(do_right)then
call htilde_mu_mat_opt_bi_ortho_tot(tmp_det,tmp_det2,$N_int,hij)
else
call htilde_mu_mat_opt_bi_ortho_tot(tmp_det2,tmp_det,$N_int,hij)
endif
!DIR$ LOOP COUNT AVG(4)
do l=1,N_st
v_t(l,k_a) = v_t(l,k_a) + hij * utl(l,kk+1)
enddo
enddo
enddo
! Diagonal contribution
! =====================
! Initial determinant is at k_a in alpha-major representation
! -----------------------------------------------------------------------
if (u_is_sparse) then
umax = 0.d0
do l=1,N_st
umax = max(umax, dabs(u_t(l,k_a)))
enddo
else
umax = 1.d0
endif
if (umax < 1.d-20) cycle
krow = psi_bilinear_matrix_rows(k_a)
ASSERT (krow <= N_det_alpha_unique)
kcol = psi_bilinear_matrix_columns(k_a)
ASSERT (kcol <= N_det_beta_unique)
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)
double precision, external :: diag_H_mat_elem
double precision :: hmono, htwoe, hthree
! hij = diag_H_mat_elem(tmp_det,$N_int)
call diag_htilde_mu_mat_fock_bi_ortho ($N_int, tmp_det, hmono, htwoe, hthree, hij)
call get_s2(tmp_det,tmp_det,$N_int,sij)
!DIR$ LOOP COUNT AVG(4)
do l=1,N_st
v_t(l,k_a) = v_t(l,k_a) + hij * u_t(l,k_a)
s_t(l,k_a) = s_t(l,k_a) + sij * u_t(l,k_a)
enddo
end do
!$OMP END DO
deallocate(buffer, singles_a, singles_b, doubles, idx, utl)
!$OMP END PARALLEL
end
SUBST [ N_int ]
1;;
2;;
3;;
4;;
N_int;;
END_TEMPLATE

3
src/tc_bi_ortho/u0_h_u0.irp.f → src/tc_bi_ortho/h_tc_u0.irp.f
View File

@ -93,9 +93,6 @@ subroutine H_tc_u_0_nstates_openmp(v_0,u_0,N_st,sze, do_right)
double precision, allocatable :: u_t(:,:), v_t(:,:) double precision, allocatable :: u_t(:,:), v_t(:,:)
!DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: u_t !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: u_t
allocate(u_t(N_st,N_det),v_t(N_st,N_det)) allocate(u_t(N_st,N_det),v_t(N_st,N_det))
! provide mo_bi_ortho_tc_one_e mo_bi_ortho_tc_two_e
! provide ref_tc_energy_tot fock_op_2_e_tc_closed_shell
! provide eff_2_e_from_3_e_ab eff_2_e_from_3_e_aa eff_2_e_from_3_e_bb
do k=1,N_st do k=1,N_st
call dset_order(u_0(1,k),psi_bilinear_matrix_order,N_det) call dset_order(u_0(1,k),psi_bilinear_matrix_order,N_det)
enddo enddo

3
src/tc_bi_ortho/tc_bi_ortho.irp.f
View File

@ -7,9 +7,6 @@ program tc_bi_ortho
! !
END_DOC END_DOC
implicit none
print *, 'Hello world'
my_grid_becke = .True. my_grid_becke = .True.
my_n_pt_r_grid = 30 my_n_pt_r_grid = 30
my_n_pt_a_grid = 50 my_n_pt_a_grid = 50

220
src/tc_bi_ortho/tc_h_eigvectors.irp.f
View File

@ -25,8 +25,6 @@ subroutine diagonalize_CI_tc
psi_r_coef_bi_ortho(i,j) = reigvec_tc_bi_orth(i,j) psi_r_coef_bi_ortho(i,j) = reigvec_tc_bi_orth(i,j)
enddo enddo
enddo enddo
! psi_energy(1:N_states) = CI_electronic_energy(1:N_states)
! psi_s2(1:N_states) = CI_s2(1:N_states)
SOFT_TOUCH psi_l_coef_bi_ortho psi_r_coef_bi_ortho SOFT_TOUCH psi_l_coef_bi_ortho psi_r_coef_bi_ortho
end end
@ -37,6 +35,7 @@ end
&BEGIN_PROVIDER [double precision, eigval_left_tc_bi_orth, (N_states)] &BEGIN_PROVIDER [double precision, eigval_left_tc_bi_orth, (N_states)]
&BEGIN_PROVIDER [double precision, reigvec_tc_bi_orth, (N_det,N_states)] &BEGIN_PROVIDER [double precision, reigvec_tc_bi_orth, (N_det,N_states)]
&BEGIN_PROVIDER [double precision, leigvec_tc_bi_orth, (N_det,N_states)] &BEGIN_PROVIDER [double precision, leigvec_tc_bi_orth, (N_det,N_states)]
&BEGIN_PROVIDER [double precision, s2_eigvec_tc_bi_orth, (N_states)]
&BEGIN_PROVIDER [double precision, norm_ground_left_right_bi_orth ] &BEGIN_PROVIDER [double precision, norm_ground_left_right_bi_orth ]
BEGIN_DOC BEGIN_DOC
@ -47,76 +46,163 @@ end
integer :: i, idx_dress, j, istate, k integer :: i, idx_dress, j, istate, k
logical :: converged, dagger logical :: converged, dagger
integer :: n_real_tc_bi_orth_eigval_right,igood_r,igood_l integer :: n_real_tc_bi_orth_eigval_right,igood_r,igood_l
integer, allocatable :: iorder(:) double precision, allocatable :: reigvec_tc_bi_orth_tmp(:,:),leigvec_tc_bi_orth_tmp(:,:),eigval_right_tmp(:)
double precision, allocatable :: reigvec_tc_bi_orth_tmp(:,:), leigvec_tc_bi_orth_tmp(:,:), eigval_right_tmp(:) double precision, allocatable :: s2_values_tmp(:), H_prime(:,:), expect_e(:)
double precision, allocatable :: coef_hf_r(:),coef_hf_l(:), Stmp(:,:) double precision, parameter :: alpha = 0.1d0
integer :: i_good_state,i_other_state, i_state
integer, allocatable :: index_good_state_array(:)
logical, allocatable :: good_state_array(:)
double precision, allocatable :: coef_hf_r(:),coef_hf_l(:)
double precision, allocatable :: Stmp(:,:)
integer, allocatable :: iorder(:)
PROVIDE N_det N_int PROVIDE N_det N_int
if(n_det .le. N_det_max_full) then if(n_det.le.N_det_max_full)then
allocate(reigvec_tc_bi_orth_tmp(N_det,N_det),leigvec_tc_bi_orth_tmp(N_det,N_det),eigval_right_tmp(N_det),expect_e(N_det))
allocate(reigvec_tc_bi_orth_tmp(N_det,N_det), leigvec_tc_bi_orth_tmp(N_det,N_det), eigval_right_tmp(N_det)) allocate (H_prime(N_det,N_det),s2_values_tmp(N_det))
H_prime(1:N_det,1:N_det) = htilde_matrix_elmt_bi_ortho(1:N_det,1:N_det)
call non_hrmt_real_diag( N_det, htilde_matrix_elmt_bi_ortho & if(s2_eig)then
, leigvec_tc_bi_orth_tmp, reigvec_tc_bi_orth_tmp, n_real_tc_bi_orth_eigval_right, eigval_right_tmp) H_prime(1:N_det,1:N_det) += alpha * S2_matrix_all_dets(1:N_det,1:N_det)
do j=1,N_det
allocate(coef_hf_r(N_det), coef_hf_l(N_det), iorder(N_det)) H_prime(j,j) = H_prime(j,j) - alpha*expected_s2
do i = 1, N_det enddo
iorder(i) = i
coef_hf_r(i) = -dabs(reigvec_tc_bi_orth_tmp(index_HF_psi_det,i))
enddo
call dsort(coef_hf_r, iorder, N_det)
igood_r = iorder(1)
print*, 'igood_r, coef_hf_r = ', igood_r, coef_hf_r(1)
do i = 1, N_det
iorder(i) = i
coef_hf_l(i) = -dabs(leigvec_tc_bi_orth_tmp(index_HF_psi_det,i))
enddo
call dsort(coef_hf_l, iorder, N_det)
igood_l = iorder(1)
print*, 'igood_l, coef_hf_l = ', igood_l, coef_hf_l(1)
if(igood_r .ne. igood_l .and. igood_r .ne. 1)then
print *,''
print *,'Warning, the left and right eigenvectors are "not the same" '
print *,'Warning, the ground state is not dominated by HF...'
print *,'State with largest RIGHT coefficient of HF ',igood_r
print *,'coef of HF in RIGHT eigenvector = ',reigvec_tc_bi_orth_tmp(index_HF_psi_det,igood_r)
print *,'State with largest LEFT coefficient of HF ',igood_l
print *,'coef of HF in LEFT eigenvector = ',leigvec_tc_bi_orth_tmp(index_HF_psi_det,igood_l)
endif endif
call non_hrmt_real_diag(N_det,H_prime,&
leigvec_tc_bi_orth_tmp,reigvec_tc_bi_orth_tmp,&
n_real_tc_bi_orth_eigval_right,eigval_right_tmp)
! do i = 1, N_det
! call get_H_tc_s2_l0_r0(leigvec_tc_bi_orth_tmp(1,i),reigvec_tc_bi_orth_tmp(1,i),1,N_det,expect_e(i), s2_values_tmp(i))
! enddo
call get_H_tc_s2_l0_r0(leigvec_tc_bi_orth_tmp,reigvec_tc_bi_orth_tmp,N_det,N_det,expect_e, s2_values_tmp)
allocate(index_good_state_array(N_det),good_state_array(N_det))
i_state = 0
good_state_array = .False.
if(s2_eig)then
if (only_expected_s2) then
do j=1,N_det
! Select at least n_states states with S^2 values closed to "expected_s2"
! print*,'s2_values_tmp(j) = ',s2_values_tmp(j),eigval_right_tmp(j),expect_e(j)
if(dabs(s2_values_tmp(j)-expected_s2).le.0.5d0)then
i_state +=1
index_good_state_array(i_state) = j
good_state_array(j) = .True.
endif
if(i_state.eq.N_states) then
exit
endif
enddo
else
do j=1,N_det
index_good_state_array(j) = j
good_state_array(j) = .True.
enddo
endif
if(i_state .ne.0)then
! Fill the first "i_state" states that have a correct S^2 value
do j = 1, i_state
do i=1,N_det
reigvec_tc_bi_orth(i,j) = reigvec_tc_bi_orth_tmp(i,index_good_state_array(j))
leigvec_tc_bi_orth(i,j) = leigvec_tc_bi_orth_tmp(i,index_good_state_array(j))
enddo
eigval_right_tc_bi_orth(j) = expect_e(index_good_state_array(j))
eigval_left_tc_bi_orth(j) = expect_e(index_good_state_array(j))
s2_eigvec_tc_bi_orth(j) = s2_values_tmp(index_good_state_array(j))
enddo
i_other_state = 0
do j = 1, N_det
if(good_state_array(j))cycle
i_other_state +=1
if(i_state+i_other_state.gt.n_states)then
exit
endif
do i=1,N_det
reigvec_tc_bi_orth(i,i_state+i_other_state) = reigvec_tc_bi_orth_tmp(i,j)
leigvec_tc_bi_orth(i,i_state+i_other_state) = leigvec_tc_bi_orth_tmp(i,j)
enddo
eigval_right_tc_bi_orth(i_state+i_other_state) = eigval_right_tmp(j)
eigval_left_tc_bi_orth (i_state+i_other_state) = eigval_right_tmp(j)
s2_eigvec_tc_bi_orth(i_state+i_other_state) = s2_values_tmp(i_state+i_other_state)
enddo
else ! istate == 0
print*,''
print*,'!!!!!!!! WARNING !!!!!!!!!'
print*,' Within the ',N_det,'determinants selected'
print*,' and the ',N_states_diag,'states requested'
print*,' We did not find only states with S^2 values close to ',expected_s2
print*,' We will then set the first N_states eigenvectors of the H matrix'
print*,' as the CI_eigenvectors'
print*,' You should consider more states and maybe ask for s2_eig to be .True. or just enlarge the CI space'
print*,''
do j=1,min(N_states_diag,N_det)
do i=1,N_det
leigvec_tc_bi_orth(i,j) = leigvec_tc_bi_orth_tmp(i,j)
reigvec_tc_bi_orth(i,j) = reigvec_tc_bi_orth_tmp(i,j)
enddo
eigval_right_tc_bi_orth(j) = eigval_right_tmp(j)
eigval_left_tc_bi_orth (j) = eigval_right_tmp(j)
s2_eigvec_tc_bi_orth(j) = s2_values_tmp(j)
enddo
endif ! istate .ne. 0
if(state_following_tc) then else ! s2_eig
allocate(coef_hf_r(N_det),coef_hf_l(N_det),iorder(N_det))
print *,'Following the states with the largest coef on HF' do i = 1,N_det
print *,'igood_r,igood_l',igood_r,igood_l iorder(i) = i
i= igood_r coef_hf_r(i) = -dabs(reigvec_tc_bi_orth_tmp(index_HF_psi_det,i))
eigval_right_tc_bi_orth(1) = eigval_right_tmp(i) enddo
do j = 1, N_det call dsort(coef_hf_r,iorder,N_det)
reigvec_tc_bi_orth(j,1) = reigvec_tc_bi_orth_tmp(j,i) igood_r = iorder(1)
enddo print*,'igood_r, coef_hf_r = ',igood_r,coef_hf_r(1)
i= igood_l do i = 1,N_det
eigval_left_tc_bi_orth(1) = eigval_right_tmp(i) iorder(i) = i
do j = 1, N_det coef_hf_l(i) = -dabs(leigvec_tc_bi_orth_tmp(index_HF_psi_det,i))
leigvec_tc_bi_orth(j,1) = leigvec_tc_bi_orth_tmp(j,i) enddo
enddo call dsort(coef_hf_l,iorder,N_det)
igood_l = iorder(1)
else print*,'igood_l, coef_hf_l = ',igood_l,coef_hf_l(1)
do i = 1, N_states if(igood_r.ne.igood_l.and.igood_r.ne.1)then
eigval_right_tc_bi_orth(i) = eigval_right_tmp(i) print *,''
eigval_left_tc_bi_orth(i) = eigval_right_tmp(i) print *,'Warning, the left and right eigenvectors are "not the same" '
print *,'Warning, the ground state is not dominated by HF...'
print *,'State with largest RIGHT coefficient of HF ',igood_r
print *,'coef of HF in RIGHT eigenvector = ',reigvec_tc_bi_orth_tmp(index_HF_psi_det,igood_r)
print *,'State with largest LEFT coefficient of HF ',igood_l
print *,'coef of HF in LEFT eigenvector = ',leigvec_tc_bi_orth_tmp(index_HF_psi_det,igood_l)
endif
if(state_following_tc)then
print *,'Following the states with the largest coef on HF'
print *,'igood_r,igood_l',igood_r,igood_l
i= igood_r
eigval_right_tc_bi_orth(1) = eigval_right_tmp(i)
do j = 1, N_det do j = 1, N_det
reigvec_tc_bi_orth(j,i) = reigvec_tc_bi_orth_tmp(j,i) reigvec_tc_bi_orth(j,1) = reigvec_tc_bi_orth_tmp(j,i)
leigvec_tc_bi_orth(j,i) = leigvec_tc_bi_orth_tmp(j,i) ! print*,reigvec_tc_bi_orth(j,1)
enddo enddo
enddo i= igood_l
eigval_left_tc_bi_orth(1) = eigval_right_tmp(i)
do j = 1, N_det
leigvec_tc_bi_orth(j,1) = leigvec_tc_bi_orth_tmp(j,i)
enddo
else
do i = 1, N_states
eigval_right_tc_bi_orth(i) = eigval_right_tmp(i)
eigval_left_tc_bi_orth(i) = eigval_right_tmp(i)
do j = 1, N_det
reigvec_tc_bi_orth(j,i) = reigvec_tc_bi_orth_tmp(j,i)
leigvec_tc_bi_orth(j,i) = leigvec_tc_bi_orth_tmp(j,i)
enddo
enddo
endif
! check bi-orthogonality ! check bi-orthogonality
allocate(Stmp(N_states,N_states)) allocate(Stmp(N_states,N_states))
call dgemm( 'T', 'N', N_states, N_states, N_det, 1.d0 & call dgemm( 'T', 'N', N_states, N_states, N_det, 1.d0 &
, leigvec_tc_bi_orth(1,1), size(leigvec_tc_bi_orth, 1), reigvec_tc_bi_orth(1,1), size(reigvec_tc_bi_orth, 1) & , leigvec_tc_bi_orth(1,1), size(leigvec_tc_bi_orth, 1), reigvec_tc_bi_orth(1,1), size(reigvec_tc_bi_orth, 1) &
, 0.d0, Stmp, size(Stmp, 1) ) , 0.d0, Stmp(1,1), size(Stmp, 1) )
print *, ' overlap matrix between states:' print *, ' overlap matrix between states:'
do i = 1, N_states do i = 1, N_states
write(*,'(1000(F16.10,X))') Stmp(i,:) write(*,'(1000(F16.10,X))') Stmp(i,:)
@ -132,6 +218,8 @@ end
external htcdag_bi_ortho_calc_tdav external htcdag_bi_ortho_calc_tdav
external H_tc_u_0_opt external H_tc_u_0_opt
external H_tc_dagger_u_0_opt external H_tc_dagger_u_0_opt
external H_tc_s2_dagger_u_0_opt
external H_tc_s2_u_0_opt
allocate(H_jj(N_det),vec_tmp(N_det,n_states_diag)) allocate(H_jj(N_det),vec_tmp(N_det,n_states_diag))
do i = 1, N_det do i = 1, N_det
call htilde_mu_mat_bi_ortho_tot(psi_det(1,1,i), psi_det(1,1,i), N_int, H_jj(i)) call htilde_mu_mat_bi_ortho_tot(psi_det(1,1,i), psi_det(1,1,i), N_int, H_jj(i))
@ -146,7 +234,8 @@ end
vec_tmp(istate,istate) = 1.d0 vec_tmp(istate,istate) = 1.d0
enddo enddo
! call davidson_general_ext_rout_nonsym_b1space(vec_tmp, H_jj, eigval_left_tc_bi_orth, N_det, n_states, n_states_diag, converged, htcdag_bi_ortho_calc_tdav) ! call davidson_general_ext_rout_nonsym_b1space(vec_tmp, H_jj, eigval_left_tc_bi_orth, N_det, n_states, n_states_diag, converged, htcdag_bi_ortho_calc_tdav)
call davidson_general_ext_rout_nonsym_b1space(vec_tmp, H_jj, eigval_left_tc_bi_orth, N_det, n_states, n_states_diag, converged, H_tc_dagger_u_0_opt) ! call davidson_general_ext_rout_nonsym_b1space(vec_tmp, H_jj, eigval_left_tc_bi_orth, N_det, n_states, n_states_diag, converged, H_tc_dagger_u_0_opt)
call davidson_hs2_nonsym_b1space(vec_tmp, H_jj, s2_eigvec_tc_bi_orth, eigval_left_tc_bi_orth, N_det, n_states, n_states_diag, converged, H_tc_s2_dagger_u_0_opt)
do istate = 1, N_states do istate = 1, N_states
leigvec_tc_bi_orth(1:N_det,istate) = vec_tmp(1:N_det,istate) leigvec_tc_bi_orth(1:N_det,istate) = vec_tmp(1:N_det,istate)
enddo enddo
@ -161,7 +250,8 @@ end
vec_tmp(istate,istate) = 1.d0 vec_tmp(istate,istate) = 1.d0
enddo enddo
! call davidson_general_ext_rout_nonsym_b1space(vec_tmp, H_jj, eigval_right_tc_bi_orth, N_det, n_states, n_states_diag, converged, htc_bi_ortho_calc_tdav) ! call davidson_general_ext_rout_nonsym_b1space(vec_tmp, H_jj, eigval_right_tc_bi_orth, N_det, n_states, n_states_diag, converged, htc_bi_ortho_calc_tdav)
call davidson_general_ext_rout_nonsym_b1space(vec_tmp, H_jj, eigval_right_tc_bi_orth, N_det, n_states, n_states_diag, converged, H_tc_u_0_opt) ! call davidson_general_ext_rout_nonsym_b1space(vec_tmp, H_jj, eigval_right_tc_bi_orth, N_det, n_states, n_states_diag, converged, H_tc_u_0_opt)
call davidson_hs2_nonsym_b1space(vec_tmp, H_jj, s2_eigvec_tc_bi_orth, eigval_right_tc_bi_orth, N_det, n_states, n_states_diag, converged, H_tc_s2_dagger_u_0_opt)
do istate = 1, N_states do istate = 1, N_states
reigvec_tc_bi_orth(1:N_det,istate) = vec_tmp(1:N_det,istate) reigvec_tc_bi_orth(1:N_det,istate) = vec_tmp(1:N_det,istate)
enddo enddo
@ -176,7 +266,9 @@ end
do j = 1, N_det do j = 1, N_det
norm_ground_left_right_bi_orth += leigvec_tc_bi_orth(j,i) * reigvec_tc_bi_orth(j,i) norm_ground_left_right_bi_orth += leigvec_tc_bi_orth(j,i) * reigvec_tc_bi_orth(j,i)
enddo enddo
print*,'norm l/r = ',norm_ground_left_right_bi_orth print*,' state ', i
print*,' norm l/r = ', norm_ground_left_right_bi_orth
print*,' <S2> = ', s2_eigvec_tc_bi_orth(i)
enddo enddo
! --- ! ---
@ -200,8 +292,6 @@ end
deallocate(buffer) deallocate(buffer)
! ---
END_PROVIDER END_PROVIDER

159
src/tc_bi_ortho/test_s2_tc.irp.f Normal file
View File

@ -0,0 +1,159 @@
program test_tc
implicit none
read_wf = .True.
my_grid_becke = .True.
my_n_pt_r_grid = 30
my_n_pt_a_grid = 50
read_wf = .True.
touch read_wf
touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid
call routine_test_s2
call routine_test_s2_davidson
end
subroutine routine_test_s2
implicit none
logical :: do_right
integer :: sze ,i, N_st, j
double precision :: sij, accu_e, accu_s, accu_e_0, accu_s_0
double precision, allocatable :: v_0_ref(:,:),u_0(:,:),s_0_ref(:,:)
double precision, allocatable :: v_0_new(:,:),s_0_new(:,:)
sze = N_det
N_st = 1
allocate(v_0_ref(N_det,1),u_0(N_det,1),s_0_ref(N_det,1),s_0_new(N_det,1),v_0_new(N_det,1))
print*,'Checking first the Left '
do_right = .False.
do i = 1, sze
u_0(i,1) = psi_l_coef_bi_ortho(i,1)
enddo
call H_tc_u_0_nstates_openmp(v_0_ref,u_0,N_st,sze, do_right)
s_0_ref = 0.d0
do i = 1, sze
do j = 1, sze
call get_s2(psi_det(1,1,i),psi_det(1,1,j),N_int,sij)
s_0_ref(i,1) += u_0(j,1) * sij
enddo
enddo
call H_tc_s2_u_0_nstates_openmp(v_0_new,s_0_new,u_0,N_st,sze, do_right)
accu_e = 0.d0
accu_s = 0.d0
accu_e_0 = 0.d0
accu_s_0 = 0.d0
do i = 1, sze
accu_e_0 += v_0_ref(i,1) * psi_r_coef_bi_ortho(i,1)
accu_s_0 += s_0_ref(i,1) * psi_r_coef_bi_ortho(i,1)
accu_e += dabs(v_0_ref(i,1) - v_0_new(i,1))
accu_s += dabs(s_0_ref(i,1) - s_0_new(i,1))
enddo
print*,'accu_e = ',accu_e
print*,'accu_s = ',accu_s
print*,'accu_e_0 = ',accu_e_0
print*,'accu_s_0 = ',accu_s_0
print*,'Checking then the right '
do_right = .True.
do i = 1, sze
u_0(i,1) = psi_r_coef_bi_ortho(i,1)
enddo
call H_tc_u_0_nstates_openmp(v_0_ref,u_0,N_st,sze, do_right)
s_0_ref = 0.d0
do i = 1, sze
do j = 1, sze
call get_s2(psi_det(1,1,i),psi_det(1,1,j),N_int,sij)
s_0_ref(i,1) += u_0(j,1) * sij
enddo
enddo
call H_tc_s2_u_0_nstates_openmp(v_0_new,s_0_new,u_0,N_st,sze, do_right)
accu_e = 0.d0
accu_s = 0.d0
accu_e_0 = 0.d0
accu_s_0 = 0.d0
do i = 1, sze
accu_e_0 += v_0_ref(i,1) * psi_l_coef_bi_ortho(i,1)
accu_s_0 += s_0_ref(i,1) * psi_l_coef_bi_ortho(i,1)
accu_e += dabs(v_0_ref(i,1) - v_0_new(i,1))
accu_s += dabs(s_0_ref(i,1) - s_0_new(i,1))
enddo
print*,'accu_e = ',accu_e
print*,'accu_s = ',accu_s
print*,'accu_e_0 = ',accu_e_0
print*,'accu_s_0 = ',accu_s_0
end
subroutine routine_test_s2_davidson
implicit none
double precision, allocatable :: H_jj(:),vec_tmp(:,:), energies(:) , s2(:)
integer :: i,istate
logical :: converged
external H_tc_s2_dagger_u_0_opt
external H_tc_s2_u_0_opt
allocate(H_jj(N_det),vec_tmp(N_det,n_states_diag),energies(n_states_diag), s2(n_states_diag))
do i = 1, N_det
call htilde_mu_mat_bi_ortho_tot(psi_det(1,1,i), psi_det(1,1,i), N_int, H_jj(i))
enddo
! Preparing the left-eigenvector
print*,'Computing the left-eigenvector '
vec_tmp = 0.d0
do istate = 1, N_states
vec_tmp(1:N_det,istate) = psi_l_coef_bi_ortho(1:N_det,istate)
enddo
do istate = N_states+1, n_states_diag
vec_tmp(istate,istate) = 1.d0
enddo
do istate = 1, N_states
leigvec_tc_bi_orth(1:N_det,istate) = vec_tmp(1:N_det,istate)
enddo
call davidson_hs2_nonsym_b1space(vec_tmp, H_jj, s2, energies, N_det, n_states, n_states_diag, converged, H_tc_s2_dagger_u_0_opt)
double precision, allocatable :: v_0_new(:,:),s_0_new(:,:)
integer :: sze,N_st
logical :: do_right
sze = N_det
N_st = 1
do_right = .False.
allocate(s_0_new(N_det,1),v_0_new(N_det,1))
call H_tc_s2_u_0_nstates_openmp(v_0_new,s_0_new,vec_tmp,N_st,sze, do_right)
double precision :: accu_e_0, accu_s_0
accu_e_0 = 0.d0
accu_s_0 = 0.d0
do i = 1, sze
accu_e_0 += v_0_new(i,1) * vec_tmp(i,1)
accu_s_0 += s_0_new(i,1) * vec_tmp(i,1)
enddo
print*,'energies = ',energies
print*,'s2 = ',s2
print*,'accu_e_0',accu_e_0
print*,'accu_s_0',accu_s_0
! Preparing the right-eigenvector
print*,'Computing the right-eigenvector '
vec_tmp = 0.d0
do istate = 1, N_states
vec_tmp(1:N_det,istate) = psi_r_coef_bi_ortho(1:N_det,istate)
enddo
do istate = N_states+1, n_states_diag
vec_tmp(istate,istate) = 1.d0
enddo
do istate = 1, N_states
leigvec_tc_bi_orth(1:N_det,istate) = vec_tmp(1:N_det,istate)
enddo
call davidson_hs2_nonsym_b1space(vec_tmp, H_jj, s2, energies, N_det, n_states, n_states_diag, converged, H_tc_s2_u_0_opt)
sze = N_det
N_st = 1
do_right = .True.
v_0_new = 0.d0
s_0_new = 0.d0
call H_tc_s2_u_0_nstates_openmp(v_0_new,s_0_new,vec_tmp,N_st,sze, do_right)
accu_e_0 = 0.d0
accu_s_0 = 0.d0
do i = 1, sze
accu_e_0 += v_0_new(i,1) * vec_tmp(i,1)
accu_s_0 += s_0_new(i,1) * vec_tmp(i,1)
enddo
print*,'energies = ',energies
print*,'s2 = ',s2
print*,'accu_e_0',accu_e_0
print*,'accu_s_0',accu_s_0
end

4
src/tc_scf/tc_scf.irp.f
View File

@ -11,8 +11,8 @@ program tc_scf
print *, ' starting ...' print *, ' starting ...'
my_grid_becke = .True. my_grid_becke = .True.
my_n_pt_r_grid = 60 my_n_pt_r_grid = 30
my_n_pt_a_grid = 110 my_n_pt_a_grid = 50
! my_n_pt_r_grid = 10 ! small grid for quick debug ! my_n_pt_r_grid = 10 ! small grid for quick debug
! my_n_pt_a_grid = 26 ! small grid for quick debug ! my_n_pt_a_grid = 26 ! small grid for quick debug
touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid