mirror of
https://github.com/QuantumPackage/qp2.git
synced 2024-11-15 10:33:50 +01:00
inactive --> virtual one-e term gradient ok
This commit is contained in:
eginer
parent
2d383d09c6
commit
8729a7ca45
78
src/casscf_tc_bi/grad_dm.irp.f
Normal file
78
src/casscf_tc_bi/grad_dm.irp.f
Normal file
@ -0,0 +1,78 @@
|
||||
BEGIN_PROVIDER [real*8, gradvec_tc_r, (0:3,nMonoEx)]
|
||||
&BEGIN_PROVIDER [real*8, gradvec_tc_l, (0:3,nMonoEx)]
|
||||
implicit none
|
||||
integer :: ii,tt,aa,indx
|
||||
integer :: i,t,a,fff
|
||||
double precision :: res_l(0:3), res_r(0:3)
|
||||
gradvec_tc_l = 0.d0
|
||||
gradvec_tc_r = 0.d0
|
||||
do i=1,n_core_inact_orb
|
||||
ii=list_core_inact(i)
|
||||
do t=1,n_act_orb
|
||||
tt=list_act(t)
|
||||
indx = mat_idx_c_a(i,t)
|
||||
call gradvec_tc_it(ii,tt,res_l)
|
||||
call gradvec_tc_it(tt,ii,res_r)
|
||||
do fff = 0,3
|
||||
gradvec_tc_l(fff,indx)=res_l(fff)
|
||||
gradvec_tc_r(fff,indx)=res_r(fff)
|
||||
enddo
|
||||
end do
|
||||
end do
|
||||
|
||||
do i=1,n_core_inact_orb
|
||||
ii=list_core_inact(i)
|
||||
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)
|
||||
do fff = 0,3
|
||||
gradvec_tc_l(fff,indx)=res_l(fff)
|
||||
gradvec_tc_r(fff,indx)=res_r(fff)
|
||||
enddo
|
||||
end do
|
||||
end do
|
||||
|
||||
do t=1,n_act_orb
|
||||
do a=1,n_virt_orb
|
||||
indx = mat_idx_a_v(i,a)
|
||||
! gradvec_tc_l(indx)=gradvec_ta(t,a)
|
||||
end do
|
||||
end do
|
||||
END_PROVIDER
|
||||
|
||||
subroutine gradvec_tc_ia(i,a,res)
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! doubly occupied --> virtual TC gradient
|
||||
!
|
||||
! Corresponds to <X0|H E_i^a|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)
|
||||
|
||||
end
|
||||
|
||||
subroutine gradvec_tc_it(i,t,res)
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! doubly occupied --> active TC gradient
|
||||
!
|
||||
! Corresponds to <X0|H E_i^t|Phi_0>
|
||||
END_DOC
|
||||
integer, intent(in) :: i,t
|
||||
double precision, intent(out) :: res(0:3)
|
||||
integer :: rr,r,ss,s
|
||||
double precision :: dm
|
||||
res = 0.d0
|
||||
res(1) = -2 * mo_bi_ortho_tc_one_e(i,t)
|
||||
do rr = 1, n_act_orb
|
||||
r = list_act(rr)
|
||||
dm = tc_transition_matrix_mo(t,r,1,1)
|
||||
res(1) += mo_bi_ortho_tc_one_e(i,r) * dm
|
||||
enddo
|
||||
|
||||
end
|
||||
109
src/casscf_tc_bi/grad_old.irp.f
Normal file
109
src/casscf_tc_bi/grad_old.irp.f
Normal file
@ -0,0 +1,109 @@
|
||||
|
||||
BEGIN_PROVIDER [real*8, gradvec_detail_right_old, (0:3,nMonoEx)]
|
||||
&BEGIN_PROVIDER [real*8, gradvec_detail_left_old, (0:3,nMonoEx)]
|
||||
BEGIN_DOC
|
||||
! calculate the orbital gradient <Psi| H E_pq |Psi> by hand, i.e. for
|
||||
! each determinant I we determine the string E_pq |I> (alpha and beta
|
||||
! separately) and generate <Psi|H E_pq |I>
|
||||
! sum_I c_I <Psi|H E_pq |I> is then the pq component of the orbital
|
||||
! gradient
|
||||
! E_pq = a^+_pa_q + a^+_Pa_Q
|
||||
END_DOC
|
||||
implicit none
|
||||
integer :: ii,tt,aa,indx,ihole,ipart,istate,ll
|
||||
real*8 :: res_l(0:3), res_r(0:3)
|
||||
|
||||
do ii = 1, n_core_inact_orb
|
||||
ihole = list_core_inact(ii)
|
||||
do aa = 1, n_virt_orb
|
||||
ipart = list_virt(aa)
|
||||
indx = mat_idx_c_v(ii,aa)
|
||||
call calc_grad_elem_h_tc(ihole,ipart,res_l, res_r)
|
||||
do ll = 0, 3
|
||||
gradvec_detail_left_old (ll,indx)=res_l(ll)
|
||||
gradvec_detail_right_old(ll,indx)=res_r(ll)
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
! do indx=1,nMonoEx
|
||||
! ihole=excit(1,indx)
|
||||
! ipart=excit(2,indx)
|
||||
! call calc_grad_elem_h_tc(ihole,ipart,res_l, res_r)
|
||||
! do ll = 0, 3
|
||||
! gradvec_detail_left_old (ll,indx)=res_l(ll)
|
||||
! gradvec_detail_right_old(ll,indx)=res_r(ll)
|
||||
! enddo
|
||||
! end do
|
||||
|
||||
real*8 :: norm_grad_left, norm_grad_right
|
||||
norm_grad_left=0.d0
|
||||
norm_grad_right=0.d0
|
||||
do indx=1,nMonoEx
|
||||
norm_grad_left+=gradvec_detail_left_old(0,indx)*gradvec_detail_left_old(0,indx)
|
||||
norm_grad_right+=gradvec_detail_right_old(0,indx)*gradvec_detail_right_old(0,indx)
|
||||
end do
|
||||
norm_grad_left=sqrt(norm_grad_left)
|
||||
norm_grad_right=sqrt(norm_grad_right)
|
||||
! if (bavard) then
|
||||
write(6,*)
|
||||
write(6,*) ' Norm of the LEFT orbital gradient (via <0|EH|0>) : ', norm_grad_left
|
||||
write(6,*) ' Norm of the RIGHT orbital gradient (via <0|HE|0>) : ', norm_grad_right
|
||||
write(6,*)
|
||||
! endif
|
||||
|
||||
|
||||
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_l = <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
|
||||
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
|
||||
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
|
||||
call det_extract(det_mu,mu,N_int)
|
||||
do ispin=1,2
|
||||
! do the monoexcitation on it
|
||||
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)
|
||||
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,N_det,N_states,i_H_chi_array,i_H_phi_array)
|
||||
! print*,i_H_chi_array(1,1),i_H_phi_array(1,1)
|
||||
do istate=1,N_states
|
||||
do ll = 0,3
|
||||
res_l(ll)+=i_H_chi_array(ll,istate)*psi_r_coef_bi_ortho(mu,istate)*phase
|
||||
res_r(ll)+=i_H_phi_array(ll,istate)*psi_l_coef_bi_ortho(mu,istate)*phase
|
||||
enddo
|
||||
end do
|
||||
end if
|
||||
end do
|
||||
end do
|
||||
|
||||
! state-averaged gradient
|
||||
res_l*=1.d0/dble(N_states)
|
||||
res_r*=1.d0/dble(N_states)
|
||||
|
||||
end
|
||||
|
||||
94
src/casscf_tc_bi/gradient.irp.f
Normal file
94
src/casscf_tc_bi/gradient.irp.f
Normal file
@ -0,0 +1,94 @@
|
||||
use bitmasks
|
||||
|
||||
BEGIN_PROVIDER [ integer, nMonoEx ]
|
||||
BEGIN_DOC
|
||||
! Number of single excitations
|
||||
END_DOC
|
||||
implicit none
|
||||
nMonoEx=n_core_inact_orb*n_act_orb+n_core_inact_orb*n_virt_orb+n_act_orb*n_virt_orb
|
||||
END_PROVIDER
|
||||
|
||||
BEGIN_PROVIDER [integer, n_c_a_prov]
|
||||
&BEGIN_PROVIDER [integer, n_c_v_prov]
|
||||
&BEGIN_PROVIDER [integer, n_a_v_prov]
|
||||
implicit none
|
||||
n_c_a_prov = n_core_inact_orb * n_act_orb
|
||||
n_c_v_prov = n_core_inact_orb * n_virt_orb
|
||||
n_a_v_prov = n_act_orb * n_virt_orb
|
||||
END_PROVIDER
|
||||
|
||||
BEGIN_PROVIDER [integer, excit, (2,nMonoEx)]
|
||||
&BEGIN_PROVIDER [character*3, excit_class, (nMonoEx)]
|
||||
&BEGIN_PROVIDER [integer, list_idx_c_a, (3,n_c_a_prov) ]
|
||||
&BEGIN_PROVIDER [integer, list_idx_c_v, (3,n_c_v_prov) ]
|
||||
&BEGIN_PROVIDER [integer, list_idx_a_v, (3,n_a_v_prov) ]
|
||||
&BEGIN_PROVIDER [integer, mat_idx_c_a, (n_core_inact_orb,n_act_orb)
|
||||
&BEGIN_PROVIDER [integer, mat_idx_c_v, (n_core_inact_orb,n_virt_orb)
|
||||
&BEGIN_PROVIDER [integer, mat_idx_a_v, (n_act_orb,n_virt_orb)
|
||||
BEGIN_DOC
|
||||
! a list of the orbitals involved in the excitation
|
||||
END_DOC
|
||||
|
||||
implicit none
|
||||
integer :: i,t,a,ii,tt,aa,indx,indx_tmp
|
||||
indx=0
|
||||
indx_tmp = 0
|
||||
do ii=1,n_core_inact_orb
|
||||
i=list_core_inact(ii)
|
||||
do tt=1,n_act_orb
|
||||
t=list_act(tt)
|
||||
indx+=1
|
||||
excit(1,indx)=i
|
||||
excit(2,indx)=t
|
||||
excit_class(indx)='c-a'
|
||||
indx_tmp += 1
|
||||
list_idx_c_a(1,indx_tmp) = indx
|
||||
list_idx_c_a(2,indx_tmp) = ii
|
||||
list_idx_c_a(3,indx_tmp) = tt
|
||||
mat_idx_c_a(ii,tt) = indx
|
||||
end do
|
||||
end do
|
||||
|
||||
indx_tmp = 0
|
||||
do ii=1,n_core_inact_orb
|
||||
i=list_core_inact(ii)
|
||||
do aa=1,n_virt_orb
|
||||
a=list_virt(aa)
|
||||
indx+=1
|
||||
excit(1,indx)=i
|
||||
excit(2,indx)=a
|
||||
excit_class(indx)='c-v'
|
||||
indx_tmp += 1
|
||||
list_idx_c_v(1,indx_tmp) = indx
|
||||
list_idx_c_v(2,indx_tmp) = ii
|
||||
list_idx_c_v(3,indx_tmp) = aa
|
||||
mat_idx_c_v(ii,aa) = indx
|
||||
end do
|
||||
end do
|
||||
|
||||
indx_tmp = 0
|
||||
do tt=1,n_act_orb
|
||||
t=list_act(tt)
|
||||
do aa=1,n_virt_orb
|
||||
a=list_virt(aa)
|
||||
indx+=1
|
||||
excit(1,indx)=t
|
||||
excit(2,indx)=a
|
||||
excit_class(indx)='a-v'
|
||||
indx_tmp += 1
|
||||
list_idx_a_v(1,indx_tmp) = indx
|
||||
list_idx_a_v(2,indx_tmp) = tt
|
||||
list_idx_a_v(3,indx_tmp) = aa
|
||||
mat_idx_a_v(tt,aa) = indx
|
||||
end do
|
||||
end do
|
||||
|
||||
! if (bavard) then
|
||||
write(6,*) ' Filled the table of the Monoexcitations '
|
||||
do indx=1,nMonoEx
|
||||
write(6,*) ' ex ',indx,' : ',excit(1,indx),' -> ' &
|
||||
,excit(2,indx),' ',excit_class(indx)
|
||||
end do
|
||||
! end if
|
||||
|
||||
END_PROVIDER
|
||||
Loading…
Reference in New Issue
Block a user