10
0
mirror of https://github.com/QuantumPackage/qp2.git synced 2024-12-22 20:34:58 +01:00

clarified the TC-CASSCF gradients

This commit is contained in:
eginer 2023-08-10 15:53:35 +02:00
parent 1a2632c280
commit ee2c470054
3 changed files with 45 additions and 36 deletions

View File

@ -24,8 +24,7 @@
do a=1,n_virt_orb
indx = mat_idx_c_v(i,a)
aa=list_virt(a)
call gradvec_tc_ia(ii,aa,res_l)
call gradvec_tc_ia(aa,ii,res_r)
call gradvec_tc_ia(ii,aa,res_l,res_r)
do fff = 0,3
gradvec_tc_l(fff,indx)=res_l(fff)
gradvec_tc_r(fff,indx)=res_r(fff)
@ -41,17 +40,21 @@
end do
END_PROVIDER
subroutine gradvec_tc_ia(i,a,res)
subroutine gradvec_tc_ia(i,a,res_l, res_r)
implicit none
BEGIN_DOC
! doubly occupied --> virtual TC gradient
!
! Corresponds to <X0|H E_i^a|Phi_0>
! Corresponds to res_r = <X0|H E_i^a|Phi_0>,
!
! res_l = <X0|E_a^i H|Phi_0>
END_DOC
integer, intent(in) :: i,a
double precision, intent(out) :: res(0:3)
res = 0.d0
res(1) = -2 * mo_bi_ortho_tc_one_e(i,a)
double precision, intent(out) :: res_l(0:3), res_r(0:3)
res_l = 0.d0
res_r = 0.d0
res_l(1) = -2 * mo_bi_ortho_tc_one_e(a,i)
res_r(1) = -2 * mo_bi_ortho_tc_one_e(i,a)
end

View File

@ -69,45 +69,51 @@ END_PROVIDER
subroutine calc_grad_elem_h_tc(ihole,ipart,res_l, res_r)
BEGIN_DOC
! eq 18 of Siegbahn et al, Physica Scripta 1980
! we calculate res_r = <Phi| H^tc E_pq | Psi>, and res_r = <Phi| E_qp H^tc | Psi>
! q=hole, p=particle
! res_l(0) = total matrix element
! res_l(1) = one-electron part
! res_l(2) = two-electron part
! res_l(3) = three-electron part
! Computes the gradient with respect to orbital rotation BRUT FORCE
!
! res_l = <Chi| E_qp H^tc | Phi>
!
! res_r = <Chi| H^tc E_pq | Phi>
!
! q=hole, p=particle. NOTE that on res_l it is E_qp and on res_r it is E_pq
!
! res_l(0) = total matrix element, res_l(1) = one-electron part,
!
! res_l(2) = two-electron part, res_l(3) = three-electron part
!
END_DOC
implicit none
integer, intent(in) :: ihole,ipart
double precision, intent(out) :: res_l(0:3), res_r(0:3)
integer :: mu,iii,ispin,ierr,nu,istate,ll
integer(bit_kind), allocatable :: det_mu(:,:),det_mu_ex(:,:)
real*8 :: i_H_chi_array(0:3,N_states),i_H_phi_array(0:3,N_states),phase
real*8 :: chi_H_mu_ex_array(0:3,N_states),mu_ex_H_phi_array(0:3,N_states),phase
allocate(det_mu(N_int,2))
allocate(det_mu_ex(N_int,2))
res_l=0.D0
res_r=0.D0
! print*,'in i_h_psi'
! print*,ihole,ipart
do mu=1,n_det
! get the string of the determinant
! get the string of the determinant |mu>
call det_extract(det_mu,mu,N_int)
do ispin=1,2
! do the monoexcitation on it
! do the monoexcitation on it: |det_mu_ex> = a^dagger_{p,ispin} a_{q,ispin} |mu>
call det_copy(det_mu,det_mu_ex,N_int)
call do_signed_mono_excitation(det_mu,det_mu_ex,nu &
,ihole,ipart,ispin,phase,ierr)
! |det_mu_ex> = a^dagger_{p,ispin} a_{q,ispin} |mu>
if (ierr.eq.1) then
call i_H_tc_psi_phi(det_mu_ex,psi_det,psi_l_coef_bi_ortho,psi_r_coef_bi_ortho,N_int &
,N_det,psi_det_size,N_states,i_H_chi_array,i_H_phi_array)
! print*,i_H_chi_array(1,1),i_H_phi_array(1,1)
,N_det,psi_det_size,N_states,chi_H_mu_ex_array,mu_ex_H_phi_array)
! chi_H_mu_ex_array = <Chi|H E_qp |mu >
! mu_ex_H_phi_array = <mu |E_qp H |Phi>
do istate=1,N_states
do ll = 0,3
res_l(ll)+=i_H_phi_array(ll,istate)*psi_l_coef_bi_ortho(mu,istate)*phase
res_r(ll)+=i_H_chi_array(ll,istate)*psi_r_coef_bi_ortho(mu,istate)*phase
do ll = 0,3 ! loop over the body components (1e,2e,3e)
!res_l = \sum_mu c_mu^l <mu|E_qp H |Phi> = <Chi|E_qp H |Phi>
res_l(ll)+= mu_ex_H_phi_array(ll,istate)*psi_l_coef_bi_ortho(mu,istate)*phase
!res_r = \sum_mu c_mu^r <Chi|H E_qp |mu> = <Chi|H E_qp |Phi>
res_r(ll)+= chi_H_mu_ex_array(ll,istate)*psi_r_coef_bi_ortho(mu,istate)*phase
enddo
end do
end if

View File

@ -90,7 +90,7 @@ subroutine htcdag_bi_ortho_calc_tdav_slow(v, u, N_st, sze)
end
subroutine i_H_tc_psi_phi(key,keys,coef_l,coef_r,Nint,Ndet,Ndet_max,Nstate,i_H_chi_array,i_H_phi_array)
subroutine i_H_tc_psi_phi(key,keys,coef_l,coef_r,Nint,Ndet,Ndet_max,Nstate,chi_H_i_array,i_H_phi_array)
use bitmasks
implicit none
BEGIN_DOC
@ -116,7 +116,7 @@ subroutine i_H_tc_psi_phi(key,keys,coef_l,coef_r,Nint,Ndet,Ndet_max,Nstate,i_H_c
integer(bit_kind), intent(in) :: keys(Nint,2,Ndet)
integer(bit_kind), intent(in) :: key(Nint,2)
double precision, intent(in) :: coef_l(Ndet_max,Nstate),coef_r(Ndet_max,Nstate)
double precision, intent(out) :: i_H_chi_array(0:3,Nstate),i_H_phi_array(0:3,Nstate)
double precision, intent(out) :: chi_H_i_array(0:3,Nstate),i_H_phi_array(0:3,Nstate)
integer :: i, ii,j
double precision :: phase
@ -131,7 +131,7 @@ subroutine i_H_tc_psi_phi(key,keys,coef_l,coef_r,Nint,Ndet,Ndet_max,Nstate,i_H_c
ASSERT (Ndet_max >= Ndet)
allocate(idx(0:Ndet))
i_H_chi_array = 0.d0
chi_H_i_array = 0.d0
i_H_phi_array = 0.d0
call filter_connected_i_H_psi0(keys,key,Nint,Ndet,idx)
@ -142,10 +142,10 @@ subroutine i_H_tc_psi_phi(key,keys,coef_l,coef_r,Nint,Ndet,Ndet_max,Nstate,i_H_c
! computes <Chi|H_tc|i>
!DIR$ FORCEINLINE
call htilde_mu_mat_opt_bi_ortho(keys(1,1,i), key, Nint, hmono, htwoe, hthree, htot)
i_H_chi_array(0,1) = i_H_chi_array(0,1) + coef_l(i,1)*htot
i_H_chi_array(1,1) = i_H_chi_array(1,1) + coef_l(i,1)*hmono
i_H_chi_array(2,1) = i_H_chi_array(2,1) + coef_l(i,1)*htwoe
i_H_chi_array(3,1) = i_H_chi_array(3,1) + coef_l(i,1)*hthree
chi_H_i_array(0,1) = chi_H_i_array(0,1) + coef_l(i,1)*htot
chi_H_i_array(1,1) = chi_H_i_array(1,1) + coef_l(i,1)*hmono
chi_H_i_array(2,1) = chi_H_i_array(2,1) + coef_l(i,1)*htwoe
chi_H_i_array(3,1) = chi_H_i_array(3,1) + coef_l(i,1)*hthree
! computes <i|H_tc|Phi>
!DIR$ FORCEINLINE
call htilde_mu_mat_opt_bi_ortho(key,keys(1,1,i), Nint, hmono, htwoe, hthree, htot)
@ -163,10 +163,10 @@ subroutine i_H_tc_psi_phi(key,keys,coef_l,coef_r,Nint,Ndet,Ndet_max,Nstate,i_H_c
!DIR$ FORCEINLINE
call htilde_mu_mat_opt_bi_ortho(keys(1,1,i), key, Nint, hmono, htwoe, hthree, htot)
do j = 1, Nstate
i_H_chi_array(0,j) = i_H_chi_array(0,j) + coef_l(i,j)*htot
i_H_chi_array(1,j) = i_H_chi_array(1,j) + coef_l(i,j)*hmono
i_H_chi_array(2,j) = i_H_chi_array(2,j) + coef_l(i,j)*htwoe
i_H_chi_array(3,j) = i_H_chi_array(3,j) + coef_l(i,j)*hthree
chi_H_i_array(0,j) = chi_H_i_array(0,j) + coef_l(i,j)*htot
chi_H_i_array(1,j) = chi_H_i_array(1,j) + coef_l(i,j)*hmono
chi_H_i_array(2,j) = chi_H_i_array(2,j) + coef_l(i,j)*htwoe
chi_H_i_array(3,j) = chi_H_i_array(3,j) + coef_l(i,j)*hthree
enddo
! computes <i|H_tc|Phi>
!DIR$ FORCEINLINE