From ee2c470054b76bffc680f54b960daae99ed9679d Mon Sep 17 00:00:00 2001 From: eginer Date: Thu, 10 Aug 2023 15:53:35 +0200 Subject: [PATCH] clarified the TC-CASSCF gradients --- src/casscf_tc_bi/grad_dm.irp.f | 17 +++++----- src/casscf_tc_bi/grad_old.irp.f | 42 ++++++++++++++----------- src/tc_bi_ortho/h_tc_bi_ortho_psi.irp.f | 22 ++++++------- 3 files changed, 45 insertions(+), 36 deletions(-) diff --git a/src/casscf_tc_bi/grad_dm.irp.f b/src/casscf_tc_bi/grad_dm.irp.f index 7f6155ab..1618adc6 100644 --- a/src/casscf_tc_bi/grad_dm.irp.f +++ b/src/casscf_tc_bi/grad_dm.irp.f @@ -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 +! Corresponds to res_r = , +! +! res_l = 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 diff --git a/src/casscf_tc_bi/grad_old.irp.f b/src/casscf_tc_bi/grad_old.irp.f index ea6747b1..6c976d66 100644 --- a/src/casscf_tc_bi/grad_old.irp.f +++ b/src/casscf_tc_bi/grad_old.irp.f @@ -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 = , and res_r = - ! 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 = + ! + ! res_r = + ! + ! 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 = + ! mu_ex_H_phi_array = 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 = + 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 = + res_r(ll)+= chi_H_mu_ex_array(ll,istate)*psi_r_coef_bi_ortho(mu,istate)*phase enddo end do end if diff --git a/src/tc_bi_ortho/h_tc_bi_ortho_psi.irp.f b/src/tc_bi_ortho/h_tc_bi_ortho_psi.irp.f index 8524253a..e96e738e 100644 --- a/src/tc_bi_ortho/h_tc_bi_ortho_psi.irp.f +++ b/src/tc_bi_ortho/h_tc_bi_ortho_psi.irp.f @@ -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 !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 !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 !DIR$ FORCEINLINE