mirror of
https://github.com/LCPQ/quantum_package
synced 2024-12-26 14:23:43 +01:00
217 lines
7.6 KiB
Fortran
217 lines
7.6 KiB
Fortran
subroutine CISD_SC2(dets_in,u_in,energies,diag_H_elements,dim_in,sze,N_st,Nint,convergence)
|
|
use bitmasks
|
|
implicit none
|
|
BEGIN_DOC
|
|
! CISD+SC2 method :: take off all the disconnected terms of a CISD (selected or not)
|
|
!
|
|
! dets_in : bitmasks corresponding to determinants
|
|
!
|
|
! u_in : guess coefficients on the various states. Overwritten
|
|
! on exit
|
|
!
|
|
! dim_in : leftmost dimension of u_in
|
|
!
|
|
! sze : Number of determinants
|
|
!
|
|
! N_st : Number of eigenstates
|
|
!
|
|
! Initial guess vectors are not necessarily orthonormal
|
|
END_DOC
|
|
integer, intent(in) :: dim_in, sze, N_st, Nint
|
|
integer(bit_kind), intent(in) :: dets_in(Nint,2,sze)
|
|
double precision, intent(inout) :: u_in(dim_in,N_st)
|
|
double precision, intent(out) :: energies(N_st)
|
|
double precision, intent(out) :: diag_H_elements(dim_in)
|
|
double precision, intent(in) :: convergence
|
|
ASSERT (N_st > 0)
|
|
ASSERT (sze > 0)
|
|
ASSERT (Nint > 0)
|
|
ASSERT (Nint == N_int)
|
|
integer :: iter
|
|
integer :: i,j,k,l,m
|
|
logical :: converged
|
|
double precision :: overlap(N_st,N_st)
|
|
double precision :: u_dot_v, u_dot_u
|
|
|
|
integer :: degree,N_double,index_hf
|
|
double precision :: hij_elec, e_corr_double,e_corr,diag_h_mat_elem,inv_c0
|
|
double precision :: e_corr_double_before,accu,cpu_2,cpu_1
|
|
integer,allocatable :: degree_exc(:), index_double(:)
|
|
integer :: i_ok
|
|
double precision,allocatable :: e_corr_array(:),H_jj_ref(:),H_jj_dressed(:),hij_double(:)
|
|
integer(bit_kind), allocatable :: doubles(:,:,:)
|
|
|
|
|
|
allocate (doubles(Nint,2,sze),e_corr_array(sze),H_jj_ref(sze),H_jj_dressed(sze),&
|
|
index_double(sze), degree_exc(sze), hij_double(sze))
|
|
call write_time(6)
|
|
write(6,'(A)') ''
|
|
write(6,'(A)') 'CISD SC2'
|
|
write(6,'(A)') '========'
|
|
!$OMP PARALLEL DEFAULT(NONE) &
|
|
!$OMP SHARED(sze,N_st, &
|
|
!$OMP H_jj_ref,Nint,dets_in,u_in) &
|
|
!$OMP PRIVATE(i)
|
|
|
|
!$OMP DO SCHEDULE(guided)
|
|
do i=1,sze
|
|
H_jj_ref(i) = diag_h_mat_elem(dets_in(1,1,i),Nint)
|
|
enddo
|
|
!$OMP END DO NOWAIT
|
|
!$OMP END PARALLEL
|
|
|
|
N_double = 0
|
|
e_corr = 0.d0
|
|
e_corr_double = 0.d0
|
|
do i = 1, sze
|
|
call get_excitation_degree(ref_bitmask,dets_in(1,1,i),degree,Nint)
|
|
degree_exc(i) = degree+1
|
|
if(degree==0)then
|
|
index_hf=i
|
|
else if (degree == 2)then
|
|
N_double += 1
|
|
index_double(N_double) = i
|
|
doubles(:,:,N_double) = dets_in(:,:,i)
|
|
call i_H_j(ref_bitmask,dets_in(1,1,i),Nint,hij_elec)
|
|
hij_double(N_double) = hij_elec
|
|
e_corr_array(N_double) = u_in(i,1)* hij_elec
|
|
e_corr_double += e_corr_array(N_double)
|
|
e_corr += e_corr_array(N_double)
|
|
else if (degree == 1)then
|
|
call i_H_j(ref_bitmask,dets_in(1,1,i),Nint,hij_elec)
|
|
e_corr += u_in(i,1)* hij_elec
|
|
endif
|
|
enddo
|
|
inv_c0 = 1.d0/u_in(index_hf,1)
|
|
do i = 1, N_double
|
|
e_corr_array(i) = e_corr_array(i) * inv_c0
|
|
enddo
|
|
e_corr = e_corr * inv_c0
|
|
e_corr_double = e_corr_double * inv_c0
|
|
converged = .False.
|
|
e_corr_double_before = e_corr_double
|
|
iter = 0
|
|
do while (.not.converged)
|
|
iter +=1
|
|
!$OMP PARALLEL DEFAULT(NONE) &
|
|
!$OMP PRIVATE(i,j,degree,accu) &
|
|
!$OMP SHARED(H_jj_dressed,sze,H_jj_ref,index_hf,N_int,N_double,&
|
|
!$OMP dets_in,doubles,degree_exc,e_corr_array,e_corr_double)
|
|
!$OMP DO SCHEDULE(STATIC)
|
|
do i=1,sze
|
|
H_jj_dressed(i) = H_jj_ref(i)
|
|
if (i==index_hf)cycle
|
|
accu = -e_corr_double
|
|
select case (N_int)
|
|
case (1)
|
|
do j=1,N_double
|
|
degree = &
|
|
popcnt(xor( dets_in(1,1,i),doubles(1,1,j))) + &
|
|
popcnt(xor( dets_in(1,2,i),doubles(1,2,j)))
|
|
|
|
if (degree<=ishft(degree_exc(i),1)) then
|
|
accu += e_corr_array(j)
|
|
endif
|
|
enddo
|
|
case (2)
|
|
do j=1,N_double
|
|
degree = &
|
|
popcnt(xor( dets_in(1,1,i),doubles(1,1,j))) + &
|
|
popcnt(xor( dets_in(1,2,i),doubles(1,2,j))) + &
|
|
popcnt(xor( dets_in(2,1,i),doubles(2,1,j))) + &
|
|
popcnt(xor( dets_in(2,2,i),doubles(2,2,j)))
|
|
|
|
if (degree<=ishft(degree_exc(i),1)) then
|
|
accu += e_corr_array(j)
|
|
endif
|
|
enddo
|
|
case (3)
|
|
do j=1,N_double
|
|
degree = &
|
|
popcnt(xor( dets_in(1,1,i),doubles(1,1,j))) + &
|
|
popcnt(xor( dets_in(1,2,i),doubles(1,2,j))) + &
|
|
popcnt(xor( dets_in(2,1,i),doubles(2,1,j))) + &
|
|
popcnt(xor( dets_in(2,2,i),doubles(2,2,j))) + &
|
|
popcnt(xor( dets_in(3,1,i),doubles(3,1,j))) + &
|
|
popcnt(xor( dets_in(3,2,i),doubles(3,2,j)))
|
|
|
|
if (degree<=ishft(degree_exc(i),1)) then
|
|
accu += e_corr_array(j)
|
|
endif
|
|
enddo
|
|
case default
|
|
do j=1,N_double
|
|
call get_excitation_degree(dets_in(1,1,i),doubles(1,1,j),degree,N_int)
|
|
if (degree<=degree_exc(i)) then
|
|
accu += e_corr_array(j)
|
|
endif
|
|
enddo
|
|
end select
|
|
H_jj_dressed(i) -= accu
|
|
enddo
|
|
!$OMP END DO
|
|
!$OMP END PARALLEL
|
|
|
|
if(sze<=N_det_max_jacobi)then
|
|
double precision, allocatable :: eigenvectors(:,:), eigenvalues(:),H_matrix_tmp(:,:)
|
|
allocate (H_matrix_tmp(size(H_matrix_all_dets,1),sze),eigenvalues(sze),eigenvectors(size(H_matrix_all_dets,1),sze))
|
|
do j=1,sze
|
|
do i=1,sze
|
|
H_matrix_tmp(i,j) = H_matrix_all_dets(i,j)
|
|
enddo
|
|
enddo
|
|
do i = 1,sze
|
|
H_matrix_tmp(i,i) = H_jj_dressed(i)
|
|
enddo
|
|
call lapack_diag(eigenvalues,eigenvectors, &
|
|
H_matrix_tmp,size(H_matrix_all_dets,1),sze)
|
|
do j=1,min(N_states_diag,sze)
|
|
do i=1,sze
|
|
u_in(i,j) = eigenvectors(i,j)
|
|
enddo
|
|
energies(j) = eigenvalues(j)
|
|
enddo
|
|
deallocate (H_matrix_tmp, eigenvalues, eigenvectors)
|
|
else
|
|
call davidson_diag_hjj(dets_in,u_in,H_jj_dressed,energies,dim_in,sze,N_st,N_states_diag,Nint,6)
|
|
endif
|
|
|
|
e_corr_double = 0.d0
|
|
inv_c0 = 1.d0/u_in(index_hf,1)
|
|
do i = 1, N_double
|
|
e_corr_array(i) = u_in(index_double(i),1)*inv_c0 * hij_double(i)
|
|
e_corr_double += e_corr_array(i)
|
|
enddo
|
|
write(6,'(A,I3)') 'SC2 Iteration ', iter
|
|
write(6,'(A)') '------------------'
|
|
write(6,'(A)') ''
|
|
write(6,'(A)') '===== ================'
|
|
write(6,'(A)') 'State Energy '
|
|
write(6,'(A)') '===== ================'
|
|
do i=1,N_st
|
|
write(6,'(I5,1X,F16.10)') i, energies(i)+nuclear_repulsion
|
|
enddo
|
|
write(6,'(A)') '===== ================'
|
|
write(6,'(A)') ''
|
|
call write_double(6,(e_corr_double - e_corr_double_before),&
|
|
'Delta(E_corr)')
|
|
converged = dabs(e_corr_double - e_corr_double_before) < convergence
|
|
converged = converged
|
|
if (converged) then
|
|
do i = 1, dim_in
|
|
diag_H_elements(i) = H_jj_dressed(i) - H_jj_ref(i)
|
|
enddo
|
|
exit
|
|
endif
|
|
e_corr_double_before = e_corr_double
|
|
|
|
enddo
|
|
|
|
call write_time(6)
|
|
deallocate (doubles,e_corr_array,H_jj_ref,H_jj_dressed, &
|
|
index_double, degree_exc, hij_double)
|
|
|
|
end
|
|
|
|
|