mirror of
https://github.com/LCPQ/quantum_package
synced 2024-12-23 04:43:50 +01:00
Merge pull request #145 from eginer/master
Added the diagonalize_s2 option
This commit is contained in:
Emmanuel Giner
commit
0b60c76656
@ -24,7 +24,7 @@ python:
|
||||
|
||||
script:
|
||||
- ./configure --production ./config/gfortran.cfg
|
||||
- source ./quantum_package.rc ; qp_module.py install Full_CI Hartree_Fock CAS_SD MRCC_CASSD
|
||||
- source ./quantum_package.rc ; qp_module.py install Full_CI Hartree_Fock CAS_SD MRCC_CASSD All_singles
|
||||
- source ./quantum_package.rc ; ninja
|
||||
- source ./quantum_package.rc ; cd ocaml ; make ; cd -
|
||||
- source ./quantum_package.rc ; cd tests ; bats bats/qp.bats
|
||||
|
||||
@ -18,7 +18,7 @@ IRPF90_FLAGS : --ninja --align=32
|
||||
# 0 : Deactivate
|
||||
#
|
||||
[OPTION]
|
||||
MODE : OPT ; [ OPT | PROFILE | DEBUG ] : Chooses the section below
|
||||
MODE : DEBUG ; [ OPT | PROFILE | DEBUG ] : Chooses the section below
|
||||
CACHE : 1 ; Enable cache_compile.py
|
||||
OPENMP : 1 ; Append OpenMP flags
|
||||
|
||||
|
||||
@ -62,11 +62,27 @@ program full_ci
|
||||
endif
|
||||
print *, 'N_det = ', N_det
|
||||
print *, 'N_states = ', N_states
|
||||
do k = 1, N_states
|
||||
print*,'State ',k
|
||||
print *, 'PT2 = ', pt2
|
||||
print *, 'E = ', CI_energy
|
||||
print *, 'E(before)+PT2 = ', E_CI_before+pt2
|
||||
enddo
|
||||
print *, '-----'
|
||||
E_CI_before = CI_energy
|
||||
if(N_states.gt.1)then
|
||||
print*,'Variational Energy difference'
|
||||
do i = 2, N_states
|
||||
print*,'Delta E = ',CI_energy(i) - CI_energy(1)
|
||||
enddo
|
||||
endif
|
||||
if(N_states.gt.1)then
|
||||
print*,'Variational + perturbative Energy difference'
|
||||
do i = 2, N_states
|
||||
print*,'Delta E = ',E_CI_before(i)+ pt2(i) - (E_CI_before(1) + pt2(1))
|
||||
enddo
|
||||
endif
|
||||
E_CI_before = CI_energy
|
||||
call ezfio_set_full_ci_energy(CI_energy)
|
||||
if (abort_all) then
|
||||
exit
|
||||
|
||||
@ -40,6 +40,12 @@ doc: Force the wave function to be an eigenfunction of S^2
|
||||
interface: ezfio,provider,ocaml
|
||||
default: False
|
||||
|
||||
[diagonalize_s2]
|
||||
type: logical
|
||||
doc: Diagonalize the S^2 operator within the n_states_diag states required. Notice : the vectors are sorted by increasing S^2 values.
|
||||
interface: ezfio,provider,ocaml
|
||||
default: True
|
||||
|
||||
[threshold_davidson]
|
||||
type: Threshold
|
||||
doc: Thresholds of Davidson's algorithm
|
||||
|
||||
@ -36,17 +36,36 @@ END_PROVIDER
|
||||
BEGIN_PROVIDER [ double precision, CI_electronic_energy, (N_states_diag) ]
|
||||
&BEGIN_PROVIDER [ double precision, CI_eigenvectors, (N_det,N_states_diag) ]
|
||||
&BEGIN_PROVIDER [ double precision, CI_eigenvectors_s2, (N_states_diag) ]
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! Eigenvectors/values of the CI matrix
|
||||
END_DOC
|
||||
integer :: i,j
|
||||
implicit none
|
||||
double precision :: ovrlp,u_dot_v
|
||||
integer :: i_good_state
|
||||
integer, allocatable :: index_good_state_array(:)
|
||||
logical, allocatable :: good_state_array(:)
|
||||
double precision, allocatable :: s2_values_tmp(:)
|
||||
integer :: i_other_state
|
||||
double precision, allocatable :: eigenvectors(:,:), eigenvalues(:)
|
||||
integer :: i_state
|
||||
double precision :: s2,e_0
|
||||
integer :: i,j,k
|
||||
double precision, allocatable :: s2_eigvalues(:)
|
||||
double precision, allocatable :: e_array(:)
|
||||
integer, allocatable :: iorder(:)
|
||||
|
||||
do j=1,N_states_diag
|
||||
! Guess values for the "N_states_diag" states of the CI_eigenvectors
|
||||
do j=1,min(N_states_diag,N_det)
|
||||
do i=1,N_det
|
||||
CI_eigenvectors(i,j) = psi_coef(i,j)
|
||||
enddo
|
||||
enddo
|
||||
|
||||
do j=N_det+1,N_states_diag
|
||||
do i=1,N_det
|
||||
CI_eigenvectors(i,j) = 0.d0
|
||||
enddo
|
||||
enddo
|
||||
|
||||
if (diag_algorithm == "Davidson") then
|
||||
|
||||
@ -59,36 +78,77 @@ END_PROVIDER
|
||||
|
||||
else if (diag_algorithm == "Lapack") then
|
||||
|
||||
double precision, allocatable :: eigenvectors(:,:), eigenvalues(:)
|
||||
allocate (eigenvectors(size(H_matrix_all_dets,1),N_det))
|
||||
allocate (eigenvalues(N_det))
|
||||
call lapack_diag(eigenvalues,eigenvectors, &
|
||||
H_matrix_all_dets,size(H_matrix_all_dets,1),N_det)
|
||||
CI_electronic_energy(:) = 0.d0
|
||||
do i=1,N_det
|
||||
CI_eigenvectors(i,1) = eigenvectors(i,1)
|
||||
enddo
|
||||
integer :: i_state
|
||||
double precision :: s2
|
||||
if (s2_eig) then
|
||||
i_state = 0
|
||||
allocate (s2_eigvalues(N_det))
|
||||
allocate(index_good_state_array(N_det),good_state_array(N_det))
|
||||
good_state_array = .False.
|
||||
do j=1,N_det
|
||||
call get_s2_u0(psi_det,eigenvectors(1,j),N_det,size(eigenvectors,1),s2)
|
||||
print*,'s2 = ',s2
|
||||
s2_eigvalues(j) = s2
|
||||
! Select at least n_states states with S^2 values closed to "expected_s2"
|
||||
if(dabs(s2-expected_s2).le.0.3d0)then
|
||||
i_state += 1
|
||||
do i=1,N_det
|
||||
CI_eigenvectors(i,i_state) = eigenvectors(i,j)
|
||||
enddo
|
||||
CI_electronic_energy(i_state) = eigenvalues(j)
|
||||
CI_eigenvectors_s2(i_state) = s2
|
||||
i_state +=1
|
||||
index_good_state_array(i_state) = j
|
||||
good_state_array(j) = .True.
|
||||
endif
|
||||
if (i_state.ge.N_states_diag) then
|
||||
exit
|
||||
if(i_state.eq.N_states) then
|
||||
exit
|
||||
endif
|
||||
enddo
|
||||
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
|
||||
CI_eigenvectors(i,j) = eigenvectors(i,index_good_state_array(j))
|
||||
enddo
|
||||
CI_electronic_energy(j) = eigenvalues(index_good_state_array(j))
|
||||
CI_eigenvectors_s2(j) = s2_eigvalues(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_diag)then
|
||||
exit
|
||||
endif
|
||||
call get_s2_u0(psi_det,eigenvectors(1,j),N_det,size(eigenvectors,1),s2)
|
||||
do i=1,N_det
|
||||
CI_eigenvectors(i,i_state+i_other_state) = eigenvectors(i,j)
|
||||
enddo
|
||||
CI_electronic_energy(i_state+i_other_state) = eigenvalues(j)
|
||||
CI_eigenvectors_s2(i_state+i_other_state) = s2
|
||||
enddo
|
||||
|
||||
deallocate(index_good_state_array,good_state_array)
|
||||
|
||||
else
|
||||
print*,''
|
||||
print*,'!!!!!!!! WARNING !!!!!!!!!'
|
||||
print*,' Within the ',N_det,'determinants selected'
|
||||
print*,' and the ',N_states_diag,'states requested'
|
||||
print*,' We did not find any state 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 diagonalize_s2 to be .True. or just enlarge the CI space'
|
||||
print*,''
|
||||
do j=1,min(N_states_diag,N_det)
|
||||
do i=1,N_det
|
||||
CI_eigenvectors(i,j) = eigenvectors(i,j)
|
||||
enddo
|
||||
CI_electronic_energy(j) = eigenvalues(j)
|
||||
CI_eigenvectors_s2(j) = s2_eigvalues(j)
|
||||
enddo
|
||||
endif
|
||||
deallocate(s2_eigvalues)
|
||||
else
|
||||
do j=1,N_states_diag
|
||||
! Select the "N_states_diag" states of lowest energy
|
||||
do j=1,min(N_det,N_states_diag)
|
||||
call get_s2_u0(psi_det,eigenvectors(1,j),N_det,N_det,s2)
|
||||
do i=1,N_det
|
||||
CI_eigenvectors(i,j) = eigenvectors(i,j)
|
||||
@ -99,7 +159,100 @@ END_PROVIDER
|
||||
endif
|
||||
deallocate(eigenvectors,eigenvalues)
|
||||
endif
|
||||
|
||||
if(diagonalize_s2.and.n_states_diag > 1.and. n_det >= n_states_diag)then
|
||||
! Diagonalizing S^2 within the "n_states_diag" states found
|
||||
allocate(s2_eigvalues(N_states_diag))
|
||||
call diagonalize_s2_betweenstates(psi_det,CI_eigenvectors,n_det,size(psi_det,3),size(CI_eigenvectors,1),min(n_states_diag,n_det),s2_eigvalues)
|
||||
|
||||
do j = 1, N_states_diag
|
||||
do i = 1, N_det
|
||||
psi_coef(i,j) = CI_eigenvectors(i,j)
|
||||
enddo
|
||||
enddo
|
||||
|
||||
if(s2_eig)then
|
||||
|
||||
! Browsing the "n_states_diag" states and getting the lowest in energy "n_states" ones that have the S^2 value
|
||||
! closer to the "expected_s2" set as input
|
||||
|
||||
allocate(index_good_state_array(N_det),good_state_array(N_det))
|
||||
good_state_array = .False.
|
||||
i_state = 0
|
||||
do j = 1, N_states_diag
|
||||
if(dabs(s2_eigvalues(j)-expected_s2).le.0.3d0)then
|
||||
good_state_array(j) = .True.
|
||||
i_state +=1
|
||||
index_good_state_array(i_state) = j
|
||||
endif
|
||||
enddo
|
||||
! Sorting the i_state good states by energy
|
||||
allocate(e_array(i_state),iorder(i_state))
|
||||
do j = 1, i_state
|
||||
do i = 1, N_det
|
||||
CI_eigenvectors(i,j) = psi_coef(i,index_good_state_array(j))
|
||||
enddo
|
||||
CI_eigenvectors_s2(j) = s2_eigvalues(index_good_state_array(j))
|
||||
call u0_H_u_0(e_0,CI_eigenvectors(1,j),n_det,psi_det,N_int)
|
||||
CI_electronic_energy(j) = e_0
|
||||
e_array(j) = e_0
|
||||
iorder(j) = j
|
||||
enddo
|
||||
call dsort(e_array,iorder,i_state)
|
||||
do j = 1, i_state
|
||||
CI_electronic_energy(j) = e_array(j)
|
||||
CI_eigenvectors_s2(j) = s2_eigvalues(index_good_state_array(iorder(j)))
|
||||
do i = 1, N_det
|
||||
CI_eigenvectors(i,j) = psi_coef(i,index_good_state_array(iorder(j)))
|
||||
enddo
|
||||
! call u0_H_u_0(e_0,CI_eigenvectors(1,j),n_det,psi_det,N_int)
|
||||
! print*,'e = ',CI_electronic_energy(j)
|
||||
! print*,'<e> = ',e_0
|
||||
! call get_s2_u0(psi_det,CI_eigenvectors(1,j),N_det,size(CI_eigenvectors,1),s2)
|
||||
! print*,'s^2 = ',CI_eigenvectors_s2(j)
|
||||
! print*,'<s^2>= ',s2
|
||||
enddo
|
||||
deallocate(e_array,iorder)
|
||||
|
||||
! Then setting the other states without any specific energy order
|
||||
i_other_state = 0
|
||||
do j = 1, N_states_diag
|
||||
if(good_state_array(j))cycle
|
||||
i_other_state +=1
|
||||
do i = 1, N_det
|
||||
CI_eigenvectors(i,i_state + i_other_state) = psi_coef(i,j)
|
||||
enddo
|
||||
CI_eigenvectors_s2(i_state + i_other_state) = s2_eigvalues(j)
|
||||
call u0_H_u_0(e_0,CI_eigenvectors(1,i_state + i_other_state),n_det,psi_det,N_int)
|
||||
CI_electronic_energy(i_state + i_other_state) = e_0
|
||||
enddo
|
||||
deallocate(index_good_state_array,good_state_array)
|
||||
|
||||
|
||||
else
|
||||
|
||||
! Sorting the N_states_diag by energy, whatever the S^2 value is
|
||||
|
||||
allocate(e_array(n_states_diag),iorder(n_states_diag))
|
||||
do j = 1, N_states_diag
|
||||
call u0_H_u_0(e_0,CI_eigenvectors(1,j),n_det,psi_det,N_int)
|
||||
e_array(j) = e_0
|
||||
iorder(j) = j
|
||||
enddo
|
||||
call dsort(e_array,iorder,n_states_diag)
|
||||
do j = 1, N_states_diag
|
||||
CI_electronic_energy(j) = e_array(j)
|
||||
do i = 1, N_det
|
||||
CI_eigenvectors(i,j) = psi_coef(i,iorder(j))
|
||||
enddo
|
||||
CI_eigenvectors_s2(j) = s2_eigvalues(iorder(j))
|
||||
enddo
|
||||
deallocate(e_array,iorder)
|
||||
endif
|
||||
deallocate(s2_eigvalues)
|
||||
endif
|
||||
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
subroutine diagonalize_CI
|
||||
|
||||
@ -214,4 +214,175 @@ subroutine get_s2_u0(psi_keys_tmp,psi_coefs_tmp,n,nmax,s2)
|
||||
deallocate (shortcut, sort_idx, sorted, version)
|
||||
end
|
||||
|
||||
subroutine get_uJ_s2_uI(psi_keys_tmp,psi_coefs_tmp,n,nmax_coefs,nmax_keys,s2,nstates)
|
||||
implicit none
|
||||
use bitmasks
|
||||
integer(bit_kind), intent(in) :: psi_keys_tmp(N_int,2,nmax_keys)
|
||||
integer, intent(in) :: n,nmax_coefs,nmax_keys,nstates
|
||||
double precision, intent(in) :: psi_coefs_tmp(nmax_coefs,nstates)
|
||||
double precision, intent(out) :: s2(nstates,nstates)
|
||||
double precision :: s2_tmp,accu
|
||||
integer :: i,j,l,jj,ll,kk
|
||||
integer, allocatable :: idx(:)
|
||||
double precision, allocatable :: tmp(:,:)
|
||||
BEGIN_DOC
|
||||
! returns the matrix elements of S^2 "s2(i,j)" between the "nstates" states
|
||||
! psi_coefs_tmp(:,i) and psi_coefs_tmp(:,j)
|
||||
END_DOC
|
||||
s2 = 0.d0
|
||||
do ll = 1, nstates
|
||||
do jj = 1, nstates
|
||||
accu = 0.d0
|
||||
!$OMP PARALLEL DEFAULT(NONE) &
|
||||
!$OMP PRIVATE (i,j,kk,idx,tmp,s2_tmp) &
|
||||
!$OMP SHARED (ll,jj,psi_keys_tmp,psi_coefs_tmp,N_int,n,nstates) &
|
||||
!$OMP REDUCTION(+:accu)
|
||||
allocate(idx(0:n))
|
||||
!$OMP DO SCHEDULE(dynamic)
|
||||
do i = 1, n
|
||||
call get_s2(psi_keys_tmp(1,1,i),psi_keys_tmp(1,1,i),s2_tmp,N_int)
|
||||
accu += psi_coefs_tmp(i,ll) * s2_tmp * psi_coefs_tmp(i,jj)
|
||||
call filter_connected(psi_keys_tmp,psi_keys_tmp(1,1,i),N_int,i-1,idx)
|
||||
do kk=1,idx(0)
|
||||
j = idx(kk)
|
||||
call get_s2(psi_keys_tmp(1,1,i),psi_keys_tmp(1,1,j),s2_tmp,N_int)
|
||||
accu += psi_coefs_tmp(i,ll) * s2_tmp * psi_coefs_tmp(j,jj) + psi_coefs_tmp(i,jj) * s2_tmp * psi_coefs_tmp(j,ll)
|
||||
enddo
|
||||
enddo
|
||||
!$OMP END DO NOWAIT
|
||||
deallocate(idx)
|
||||
!$OMP BARRIER
|
||||
!$OMP END PARALLEL
|
||||
s2(ll,jj) += accu
|
||||
enddo
|
||||
enddo
|
||||
do i = 1, nstates
|
||||
do j =i+1,nstates
|
||||
accu = 0.5d0 * (s2(i,j) + s2(j,i))
|
||||
s2(i,j) = accu
|
||||
s2(j,i) = accu
|
||||
enddo
|
||||
enddo
|
||||
end
|
||||
|
||||
subroutine diagonalize_s2_betweenstates(keys_tmp,psi_coefs_inout,n,nmax_keys,nmax_coefs,nstates,s2_eigvalues)
|
||||
BEGIN_DOC
|
||||
! You enter with nstates vectors in psi_coefs_inout that may be coupled by S^2
|
||||
! The subroutine diagonalize the S^2 operator in the basis of these states.
|
||||
! The vectors that you obtain in output are no more coupled by S^2,
|
||||
! which does not necessary mean that they are eigenfunction of S^2.
|
||||
! n,nmax,nstates = number of determinants, physical dimension of the arrays and number of states
|
||||
! keys_tmp = array of integer(bit_kind) that represents the determinants
|
||||
! psi_coefs(i,j) = coeff of the ith determinant in the jth state
|
||||
! VECTORS ARE SUPPOSED TO BE ORTHONORMAL IN INPUT
|
||||
END_DOC
|
||||
implicit none
|
||||
use bitmasks
|
||||
integer, intent(in) :: n,nmax_keys,nmax_coefs,nstates
|
||||
integer(bit_kind), intent(in) :: keys_tmp(N_int,2,nmax_keys)
|
||||
double precision, intent(inout) :: psi_coefs_inout(nmax_coefs,nstates)
|
||||
|
||||
!integer, intent(in) :: ndets_real,ndets_keys,ndets_coefs,nstates
|
||||
!integer(bit_kind), intent(in) :: keys_tmp(N_int,2,ndets_keys)
|
||||
!double precision, intent(inout) :: psi_coefs_inout(ndets_coefs,nstates)
|
||||
double precision, intent(out) :: s2_eigvalues(nstates)
|
||||
|
||||
|
||||
double precision,allocatable :: s2(:,:),overlap(:,:)
|
||||
double precision, allocatable :: eigvalues(:),eigvectors(:,:)
|
||||
integer :: i,j,k
|
||||
double precision, allocatable :: psi_coefs_tmp(:,:)
|
||||
double precision :: accu,coef_contract
|
||||
double precision :: u_dot_u,u_dot_v
|
||||
|
||||
print*,''
|
||||
print*,'*********************************************************************'
|
||||
print*,'Cleaning the various vectors by diagonalization of the S^2 matrix ...'
|
||||
print*,''
|
||||
print*,'nstates = ',nstates
|
||||
allocate(s2(nstates,nstates),overlap(nstates,nstates))
|
||||
do i = 1, nstates
|
||||
overlap(i,i) = u_dot_u(psi_coefs_inout(1,i),n)
|
||||
do j = i+1, nstates
|
||||
overlap(i,j) = u_dot_v(psi_coefs_inout(1,j),psi_coefs_inout(1,i),n)
|
||||
overlap(j,i) = overlap(i,j)
|
||||
enddo
|
||||
enddo
|
||||
print*,'Overlap matrix in the basis of the states considered'
|
||||
do i = 1, nstates
|
||||
write(*,'(10(F16.10,X))')overlap(i,:)
|
||||
enddo
|
||||
call ortho_lowdin(overlap,size(overlap,1),nstates,psi_coefs_inout,size(psi_coefs_inout,1),n)
|
||||
print*,'passed ortho'
|
||||
|
||||
do i = 1, nstates
|
||||
overlap(i,i) = u_dot_u(psi_coefs_inout(1,i),n)
|
||||
do j = i+1, nstates
|
||||
overlap(i,j) = u_dot_v(psi_coefs_inout(1,j),psi_coefs_inout(1,i),n)
|
||||
overlap(j,i) = overlap(i,j)
|
||||
enddo
|
||||
enddo
|
||||
print*,'Overlap matrix in the basis of the Lowdin orthonormalized states '
|
||||
do i = 1, nstates
|
||||
write(*,'(10(F16.10,X))')overlap(i,:)
|
||||
enddo
|
||||
|
||||
call get_uJ_s2_uI(keys_tmp,psi_coefs_inout,n_det,size(psi_coefs_inout,1),size(keys_tmp,3),s2,nstates)
|
||||
print*,'S^2 matrix in the basis of the states considered'
|
||||
double precision :: accu_precision_diag,accu_precision_of_diag
|
||||
accu_precision_diag = 0.d0
|
||||
accu_precision_of_diag = 0.d0
|
||||
do i = 1, nstates
|
||||
do j = i+1, nstates
|
||||
if( ( dabs(s2(i,i) - s2(j,j)) .le.1.d-10 ) .and. (dabs(s2(i,j) + dabs(s2(i,j)))) .le.1.d-10) then
|
||||
s2(i,j) = 0.d0
|
||||
s2(j,i) = 0.d0
|
||||
endif
|
||||
enddo
|
||||
enddo
|
||||
do i = 1, nstates
|
||||
write(*,'(10(F10.6,X))')s2(i,:)
|
||||
enddo
|
||||
|
||||
print*,'Diagonalizing the S^2 matrix'
|
||||
|
||||
allocate(eigvalues(nstates),eigvectors(nstates,nstates))
|
||||
call lapack_diagd(eigvalues,eigvectors,s2,nstates,nstates)
|
||||
print*,'Eigenvalues of s^2'
|
||||
do i = 1, nstates
|
||||
print*,'s2 = ',eigvalues(i)
|
||||
s2_eigvalues(i) = eigvalues(i)
|
||||
enddo
|
||||
|
||||
print*,'Building the eigenvectors of the S^2 matrix'
|
||||
allocate(psi_coefs_tmp(nmax_coefs,nstates))
|
||||
psi_coefs_tmp = 0.d0
|
||||
do j = 1, nstates
|
||||
do k = 1, nstates
|
||||
coef_contract = eigvectors(k,j) ! <phi_k|Psi_j>
|
||||
do i = 1, n_det
|
||||
psi_coefs_tmp(i,j) += psi_coefs_inout(i,k) * coef_contract
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
do j = 1, nstates
|
||||
accu = 0.d0
|
||||
do i = 1, n_det
|
||||
accu += psi_coefs_tmp(i,j) * psi_coefs_tmp(i,j)
|
||||
enddo
|
||||
print*,'Norm of vector = ',accu
|
||||
accu = 1.d0/dsqrt(accu)
|
||||
do i = 1, n_det
|
||||
psi_coefs_inout(i,j) = psi_coefs_tmp(i,j) * accu
|
||||
enddo
|
||||
enddo
|
||||
!call get_uJ_s2_uI(keys_tmp,psi_coefs_inout,n_det,size(psi_coefs_inout,1),size(keys_tmp,3),s2,nstates)
|
||||
!print*,'S^2 matrix in the basis of the NEW states considered'
|
||||
!do i = 1, nstates
|
||||
! write(*,'(10(F16.10,X))')s2(i,:)
|
||||
!enddo
|
||||
|
||||
deallocate(s2,eigvalues,eigvectors,psi_coefs_tmp,overlap)
|
||||
|
||||
end
|
||||
|
||||
|
||||
@ -1507,6 +1507,33 @@ subroutine get_occ_from_key(key,occ,Nint)
|
||||
|
||||
end
|
||||
|
||||
subroutine u0_H_u_0(e_0,u_0,n,keys_tmp,Nint)
|
||||
use bitmasks
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! Computes e_0 = <u_0|H|u_0>/<u_0|u_0>
|
||||
!
|
||||
! n : number of determinants
|
||||
!
|
||||
END_DOC
|
||||
integer, intent(in) :: n,Nint
|
||||
double precision, intent(out) :: e_0
|
||||
double precision, intent(in) :: u_0(n)
|
||||
integer(bit_kind),intent(in) :: keys_tmp(Nint,2,n)
|
||||
|
||||
double precision :: H_jj(n)
|
||||
double precision :: v_0(n)
|
||||
double precision :: u_dot_u,u_dot_v,diag_H_mat_elem
|
||||
integer :: i,j
|
||||
do i = 1, n
|
||||
H_jj(i) = diag_H_mat_elem(keys_tmp(1,1,i),Nint)
|
||||
enddo
|
||||
|
||||
call H_u_0(v_0,u_0,H_jj,n,keys_tmp,Nint)
|
||||
e_0 = u_dot_v(v_0,u_0,n)/u_dot_u(u_0,n)
|
||||
end
|
||||
|
||||
|
||||
subroutine H_u_0(v_0,u_0,H_jj,n,keys_tmp,Nint)
|
||||
use bitmasks
|
||||
implicit none
|
||||
|
||||
Loading…
Reference in New Issue
Block a user