subroutine tdhf_mrsf_energy(infos)
use io_constants, only: iw
use oqp_tagarray_driver
use types, only: information
use strings, only: cstring, fstring
use basis_tools, only: basis_set
use messages, only: show_message, with_abort
use util, only: measure_time
use precision, only: dp
use int2_compute, only: int2_compute_t
use tdhf_mrsf_lib, only: int2_mrsf_data_t
use tdhf_lib, only: int2_td_data_t
use tdhf_lib, only: &
iatogen, mntoia, rparedms, rpaeig, rpavnorm, &
rpaechk, rpanewb, &
rpaprint, inivec
use tdhf_sf_lib, only: sfresvec, sfqvec, sfdmat, trfrmb, &
get_transition_density, get_transitions, &
get_transition_dipole, print_results
use tdhf_mrsf_lib, only: &
mrinivec, mrsfcbc, mrsfmntoia, mrsfesum, &
mrsfqroesum, get_mrsf_transitions, &
get_mrsf_transition_density
use mathlib, only: orthogonal_transform, orthogonal_transform_sym, &
unpack_matrix
use oqp_linalg
use printing, only: print_module_info
implicit none
character(len=*), parameter :: subroutine_name = "tdhf_mrsf_energy"
type(basis_set), pointer :: basis
type(information), target, intent(inout) :: infos
integer :: s_size, ok
real(kind=dp), allocatable :: scr2(:)
real(kind=dp), allocatable :: wrk1(:,:), qvec(:,:)
real(kind=dp), allocatable :: amo(:,:), wrk2(:,:)
real(kind=dp), allocatable :: squared_S(:)
real(kind=dp), allocatable :: amb(:,:), apb(:,:)
real(kind=dp), allocatable, target :: vl(:), vr(:)
real(kind=dp), pointer :: vl_p(:,:), vr_p(:,:)
real(kind=dp), allocatable :: xm(:), scr(:)
real(kind=dp), allocatable :: bvec_mo(:,:), for_trnsf_b_vec(:,:)
real(kind=dp), allocatable, dimension(:,:) :: fa, fb
real(kind=dp), allocatable, dimension(:) :: rnorm
real(kind=dp), allocatable, dimension(:,:,:,:) :: trden
integer, allocatable, dimension(:,:) :: trans
real(kind=dp), allocatable, target :: mrsf_density(:,:,:,:)
real(kind=dp), pointer :: fmrst2(:,:,:,:)
real(kind=dp), allocatable, target :: fmrq1(:,:,:)
real(kind=dp), allocatable :: dip(:,:,:), bvec_mo_tmp(:), eex(:)
integer :: nocca, nvira, noccb, nvirb
integer :: nbf, nbf2, xvec_dim
integer :: mxvec, ist, jst, iend, nvec, novec
integer :: iter, nv, iv, ivec
integer :: mxiter
logical :: tamm_dancoff
integer :: imax
integer :: ierr
logical :: converged
real(kind=dp) :: mxerr, cnvtol, scale_exch
integer :: maxvec, mrst, nstates, target_state
logical :: roref = .false.
logical :: debug_mode
type(int2_compute_t) :: int2_driver
type(int2_mrsf_data_t), target :: int2_data_st
type(int2_td_data_t), target :: int2_data_q
logical :: dft = .false.
integer :: scf_type, mol_mult
! tagarray
real(kind=dp), contiguous, pointer :: &
fock_a(:), dmat_a(:), mo_A(:,:), mo_energy_a(:), &
fock_b(:), dmat_b(:), mo_b(:,:), mo_energy_b(:), &
smat(:), ta(:), tb(:), td_t(:,:), bvec_mo_out(:,:), &
mrsf_energies(:)
character(len=*), parameter :: tags_alloc(3) = (/ character(len=80) :: &
OQP_td_bvec_mo, OQP_td_t, OQP_td_energies /)
character(len=*), parameter :: tags_required(9) = (/ character(len=80) :: &
OQP_FOCK_A, OQP_DM_A, OQP_E_MO_A, OQP_VEC_MO_A, OQP_FOCK_B, OQP_DM_B, OQP_E_MO_B, OQP_VEC_MO_B, OQP_SM /)
! Readings
! Files open
! 3. LOG: Write: Main output file
open (unit=iw, file=infos%log_filename, position="append")
!
call print_module_info('MRSF_TDHF_Energy','Computing Energy of MRSF-TDDFT')
! Load basis set
basis => infos%basis
basis%atoms => infos%atoms
! Input parameters
dft = infos%control%hamilton == 20 ! dft or hf
mrst = infos%tddft%mult
nstates = infos%tddft%nstate
target_state = infos%tddft%target_state
maxvec = infos%tddft%maxvec
cnvtol = infos%tddft%cnvtol
! infos%tddft%debug_mode = .True.
debug_mode = infos%tddft%debug_mode
mol_mult = infos%mol_prop%mult
if (mol_mult/=3) call show_message('MRSF-TDDFT are available for ROHF ref.&
&with ONLY triplet multiplicity(mult=3)',with_abort)
scf_type = infos%control%scftype
if (scf_type==3) roref = .true.
nbf = basis%nbf
nbf2 = nbf*(nbf+1)/2
s_size = (basis%nshell**2+basis%nshell)/2
tamm_dancoff = .true. ! tamm_dancoff: 0/1 means not doing/doing Tamm/Dancoff run
nocca = infos%mol_prop%nelec_a
nvira = nbf-nocca
noccb = infos%mol_prop%nelec_b
nvirb = nbf-noccb
if (mrst==1 .or. mrst==3 ) then
xvec_dim = nocca*nvirb
else if (mrst==5) then
xvec_dim = noccb*nvira
end if
if (mrst==1 ) then
nstates = min(nstates, xvec_dim-1)
mxvec = min(maxvec*nstates, xvec_dim-1, infos%control%maxit_dav*nstates)
else if (mrst==3) then
nstates = min(nstates, xvec_dim-3)
mxvec = min(maxvec*nstates, xvec_dim-3, infos%control%maxit_dav*nstates)
else if (mrst==5) then
nstates = min(nstates, xvec_dim)
mxvec = min(maxvec*nstates, xvec_dim, infos%control%maxit_dav*nstates)
end if
infos%tddft%nstate = nstates
nvec = min(max(nstates,6), mxvec)
call infos%dat%remove_records(tags_alloc)
call infos%dat%reserve_data(OQP_td_bvec_mo, TA_TYPE_REAL64, &
xvec_dim*nstates, (/xvec_dim, nstates/), comment=OQP_td_bvec_mo_comment)
call infos%dat%reserve_data(OQP_td_t, TA_TYPE_REAL64, nbf2*2, (/ nbf2, 2 /), comment=OQP_td_t_comment)
call infos%dat%reserve_data(OQP_td_energies, TA_TYPE_REAL64, nstates, comment=OQP_td_energies_comment)
call data_has_tags(infos%dat, tags_alloc, module_name, subroutine_name, WITH_ABORT)
call tagarray_get_data(infos%dat, OQP_td_bvec_mo, bvec_mo_out)
call tagarray_get_data(infos%dat, OQP_td_t, td_t)
call tagarray_get_data(infos%dat, OQP_td_energies, mrsf_energies)
call data_has_tags(infos%dat, tags_required, module_name, subroutine_name, WITH_ABORT)
call tagarray_get_data(infos%dat, OQP_SM, smat)
call tagarray_get_data(infos%dat, OQP_FOCK_A, fock_a)
call tagarray_get_data(infos%dat, OQP_FOCK_B, fock_b)
call tagarray_get_data(infos%dat, OQP_DM_A, dmat_a)
call tagarray_get_data(infos%dat, OQP_DM_B, dmat_b)
call tagarray_get_data(infos%dat, OQP_E_MO_A, mo_energy_a)
call tagarray_get_data(infos%dat, OQP_E_MO_B, mo_energy_b)
call tagarray_get_data(infos%dat, OQP_VEC_MO_A, mo_a)
call tagarray_get_data(infos%dat, OQP_VEC_MO_B, mo_b)
! Allocate temporary matrices for diagonalization
allocate (fa(nbf,nbf), &
fb(nbf,nbf), &
stat=ok)
if( ok/=0 ) call show_message('Cannot allocate memory',with_abort)
! Allocate TDDFT variables
allocate(xm(xvec_dim), &
bvec_mo(xvec_dim,mxvec), &
trden(nbf,nbf,nstates,nstates), &
wrk1(nbf,nbf), &
wrk2(nbf,nbf), &
amo(xvec_dim,mxvec), &
EEX(mxvec), &
squared_S(nstates), &
APB(mxvec,mxvec), &
AMB(mxvec,mxvec), &
VR(mxvec*mxvec), &
VL(mxvec*mxvec), &
for_trnsf_b_vec(mxvec,mxvec), & !
dip(3,nstates,nstates), &
scr(nbf2), &
bvec_mo_tmp(xvec_dim), &
scr2(mxvec*mxvec), &
qvec(xvec_dim,nstates), &
RNORM(nstates), &
source=0.0_dp,stat=ok)
allocate(trans(xvec_dim,2), &
source=0, stat=ok)
if( ok/=0 ) call show_message('Cannot allocate memory', with_abort)
ta => td_t(:,1)
tb => td_t(:,2)
if (mrst==1 .or. mrst==3 ) then
allocate(mrsf_density(nvec,7,nbf,nbf), &
source=0.0_dp, &
stat=ok)
else if( mrst==5 )then
allocate(fmrq1(nbf,nbf,nvec), &
source=0.0_dp, &
stat=ok)
end if
if( ok/=0 ) call show_message('Cannot allocate memory', with_abort)
scale_exch = 1.0_dp
if (infos%tddft%HFscale == -1.0_dp) &
infos%tddft%HFscale = infos%dft%HFscale
if (infos%dft%cam_flag) then
if (infos%tddft%cam_alpha == -1.0_dp) &
infos%tddft%cam_alpha = infos%dft%cam_alpha
infos%tddft%HFscale = infos%tddft%cam_alpha
if (infos%tddft%cam_beta == -1.0_dp) &
infos%tddft%cam_beta = infos%dft%cam_beta
if (infos%tddft%cam_mu == -1.0_dp) &
infos%tddft%cam_mu = infos%dft%cam_mu
end if
if (dft) scale_exch = infos%tddft%HFscale
! set spin-pair coupling
if (infos%tddft%spc_coco==-1.0_dp) &
infos%tddft%spc_coco = infos%tddft%HFscale
if (infos%tddft%spc_ovov==-1.0_dp) &
infos%tddft%spc_ovov = infos%tddft%HFscale
if (infos%tddft%spc_coov==-1.0_dp) &
infos%tddft%spc_coov = infos%tddft%HFscale
if(debug_mode.or..true.)then
write(*,'(/,5x,"Input parameters:")')
write(*,'(5x,"Number of states: ",1x,I0)') nstates
write(*,'(5x,"Number of single excitations: ",1x,I0)') xvec_dim
write(*,'(5x,"Number of atomic orbitals: ",1x,I0)') nbf
write(*,'(5x,"Number of electrons: ",1x,I0)') nocca+noccb
write(*,'(5x,"Number of occupied alpha orbitals:",1x,I0)') nocca
write(*,'(5x,"Number of occupied beta orbitals: ",1x,I0)') noccb
write(*,'(5x,"Number of virtual alpha orbitals: ",1x,I0)') nvira
write(*,'(5x,"Number of virtual beta orbitals: ",1x,I0)') nvirb
write(*,'(5x,"Maximum vectors: ",1x,I0)') mxvec
write(*,'(5x,"Initial vectors: ",1x,I0)') nvec
write(*, '(/7x,"Fitting parameters for MRSF-TDDFT")')
if (.not.infos%dft%cam_flag) then
write(*, '(10x,"Exact HF exchange:")')
write(*, '(5x,"Reference: |", t20, f6.3, t29, "|")') infos%dft%HFscale
write(*, '(5x,"Response: |", t20, f6.3, t29, "|")') infos%tddft%HFscale
else
write(*, '(10x,"CAM parametres:")')
write(*, '(16x,"| alpha | beta | mu |")')
write(*, '(5x,"Reference: |", t20, f6.3, t29, "|", t32, f6.3, t41, "|", t44, f6.3, t53, "|")') &
infos%dft%cam_alpha, infos%dft%cam_beta, infos%dft%cam_mu
write(*, '(5x,"Response: |", t20, f6.3, t29, "|", t32, f6.3, t41, "|", t44, f6.3, t53, "|")') &
infos%tddft%cam_alpha, infos%tddft%cam_beta, infos%tddft%cam_mu
end if
write(*, '(10x,"Spin-pair coupling parametres:")')
write(*, '(16x,"| CO-CO | OV-OV | CO-OV |")')
write(*, '(16x,"|", t20, f6.3, t29, "|", t32, f6.3, t41, "|", t44, f6.3, t53, "|")') &
infos%tddft%spc_coco, infos%tddft%spc_ovov, infos%tddft%spc_coov
end if
write(*,'(/,5x,46("="))')
if (mrst==1) write(*,'( 5X,"Davidson algorithm for Singlet response states")')
if (mrst==3) write(*,'( 5x,"Davidson algorithm for Triplet response states")')
if (mrst==5) write(*,'( 5x,"Davidson algorithm for Quintet response states")')
write(*,'(5x,46("="))')
! Initialize ERI (Electron Repulsion Integrals) calculations
call int2_driver%init(basis, infos)
call int2_driver%set_screening()
call flush(iw)
! Prepare for ROHF
if (roref) then
! Alapha
call orthogonal_transform_sym(nbf, nbf, fock_a, mo_a, nbf, scr)
call unpack_matrix(scr,fa)
! Beta
call orthogonal_transform_sym(nbf, nbf, fock_b, mo_b, nbf, scr)
call unpack_matrix(scr,fb)
end if
! Construct TD trial vector
if (mrst==1 .or. mrst==3) then
call mrinivec(infos, mo_energy_a, mo_energy_a, bvec_mo, xm, nvec)
else if (mrst==5) then
call inivec(mo_energy_a,mo_energy_a,bvec_mo,xm,noccb,nocca,nvec)
end if
ist = 1
iend = nvec
iter = 0
mxiter = infos%control%maxit_dav
ierr = 0
do iter = 1, mxiter
nv = iend-ist+1
if( mrst==1 .or. mrst==3 ) then
mrsf_density = 0.0_dp ! bo2v, bo1v, bco1, bco2, o21v, co12, ball
else if( mrst==5 ) then
fmrq1 = 0.0_dp
end if
do ivec = ist, iend
iv = ivec-ist+1
if (mrst==1 .or. mrst==3) then
call iatogen(bvec_mo(:,ivec), wrk1, nocca, noccb)
call mrsfcbc(infos, mo_a, mo_b, wrk1, mrsf_density(iv,:,:,:))
else if (mrst==5) then
call iatogen(bvec_mo(:,ivec), wrk1, noccb, nocca)
call orthogonal_transform('t', nbf, mo_a, wrk1, fmrq1(:,:,iv), wrk2)
end if
end do
if (mrst==1 .or. mrst==3) then
int2_data_st = int2_mrsf_data_t( &
d3 = mrsf_density(:iv,:,:,:), &
tamm_dancoff = tamm_dancoff, &
scale_exchange = scale_exch, &
scale_coulomb = scale_exch)
call int2_driver%run( &
int2_data_st, &
cam = dft.and.infos%dft%cam_flag, &
alpha = infos%tddft%cam_alpha, &
alpha_coulomb = infos%tddft%cam_alpha, &
beta = infos%tddft%cam_beta, &
beta_coulomb = infos%tddft%cam_beta, &
mu = infos%tddft%cam_mu)
fmrst2 => int2_data_st%f3(:,:,:,:,1) ! ado2v, ado1v, adco1, adco2, ao21v, aco12, agdlr
! Scaling factor if triplet
if (mrst==3) fmrst2(:,1:6,:,:) = -fmrst2(:,1:6,:,:)
! Spin pair coupling
if (infos%tddft%spc_coco /= infos%tddft%hfscale) &
fmrst2(:,6,:,:) = fmrst2(:,6,:,:) * infos%tddft%spc_coco / infos%tddft%hfscale
if (infos%tddft%spc_ovov /= infos%tddft%hfscale) &
fmrst2(:,5,:,:) = fmrst2(:,5,:,:) * infos%tddft%spc_ovov / infos%tddft%hfscale
if (infos%tddft%spc_coov /= infos%tddft%hfscale) &
fmrst2(:,1:4,:,:) = fmrst2(:,1:4,:,:) * infos%tddft%spc_coov / infos%tddft%hfscale
else if (mrst==5) then
int2_data_q = int2_td_data_t( &
d2=fmrq1(:,:,:iv), &
int_apb = .false., &
int_amb = .false., &
tamm_dancoff = tamm_dancoff, &
scale_exchange = scale_exch)
call int2_driver%run( &
int2_data_q, &
cam = dft.and.infos%dft%cam_flag, &
alpha = infos%tddft%cam_alpha, &
beta = infos%tddft%cam_beta,&
mu = infos%tddft%cam_mu)
end if
do ivec = ist, iend
iv = ivec-ist+1
if (mrst==1 .or. mrst==3) then
! Product (A-B)*X
call mrsfmntoia(infos, fmrst2(iv,:,:,:), amo, mo_a, mo_b, ivec)
call iatogen(bvec_mo(:,ivec), wrk1, nocca, noccb)
call mrsfesum(infos, wrk1, fa, fb, amo, ivec)
else if( mrst==5 )then
call mntoia(int2_data_q%amb(:,:,iv,1), amo(:,ivec), mo_a, mo_b, noccb, nocca)
! Z(I+,A-)
call iatogen(bvec_mo(:,ivec),wrk1,noccb,nocca)
! FB(I+,J+)*Z(J+,A-)
call dgemm('n','n',noccb,nbf,noccb, &
1.0_dp,fb,nbf, &
wrk1,nbf, &
0.0_dp,wrk2,noccb)
! Z(I+,B-)*FA(B-,A-)
call dgemm('n','n',noccb,nbf,nbf, &
1.0_dp,wrk1,nbf, &
fa,nbf, &
-1.0_dp,wrk2,noccb)
call mrsfqroesum(wrk2,amo, &
nocca,noccb,nbf,ivec)
end if
end do
vl_p(1:nvec, 1:nvec) => vl(1:nvec*nvec)
vr_p(1:nvec, 1:nvec) => vr(1:nvec*nvec)
call rparedms(bvec_mo,amo,amo,apb,amb,nvec,tamm_dancoff=.true.)
call rpaeig(eex,vl_p,vr_p,apb,amb,scr2,tamm_dancoff=.true.)
call rpavnorm(vr_p,vl_p,tamm_dancoff=.true.)
call rpaechk(eex,nvec,nstates,imax,tamm_dancoff=.true.)
for_trnsf_b_vec = vr_p
call sfresvec(qvec,bvec_mo,amo,vr_p,eex,nvec,rnorm,nstates)
call sfqvec(qvec,xm,eex,nstates)
call rpaprint(eex, rnorm, cnvtol, iter, imax, nstates, do_neg=.true.)
mxerr = maxval(rnorm)
! Check convergence
converged = mxerr<=cnvtol
if (converged) exit
! No space left for new vectors, exit
if (nvec==mxvec) ierr = 1
if (ierr/=0) exit
call rpanewb(nstates,bvec_mo,qvec,novec,nvec,ierr,tamm_dancoff=.true.)
! ierr=1 nvec over mxvec: not converged case
if (ierr/=0) exit
ist = novec+1
iend = nvec
end do
if (iter >= mxiter .and. .not. converged) ierr = -1
select case (ierr)
case (-1)
write(*,'(/,2X,"MRSF-TD-DFT energies NOT CONVERGED after ",I4," iterations"/)') mxiter
infos%mol_energy%Davidson_converged=.false.
case (0)
write(*,'(/,2X,"MRSF-TD-DFT energies converged in ",I4," iterations"/)') iter
infos%mol_energy%Davidson_converged=.true.
case (1)
write(*,'(/,2X,"..something is wrong.. nvec = mxvec")')
infos%mol_energy%Davidson_converged=.false.
case (2)
write(*,'(/,2x,"..something is wrong.. nvec > mxvec")')
write(*,'(3x,"nvec/mxvec =",I4,"/",I4)') nvec, mxvec
infos%mol_energy%Davidson_converged=.false.
case (3)
write(*,'(/,2x,"..something is wrong.. No vectors were added")')
infos%mol_energy%Davidson_converged=.false.
end select
call flush(iw)
call trfrmb(bvec_mo, for_trnsf_b_vec, nvec, nstates)
select case (mrst)
case(1)
do ist = 1, nstates
do jst = ist, nstates
call get_mrsf_transition_density(infos,trden(:,:,ist,jst), bvec_mo, ist, jst)
end do
end do
squared_S(:) = 0.0_dp
call get_mrsf_transitions(trans, nocca, noccb, nbf)
write(*,'(/,2x,35("="),/,2x,&
&"Spin-adapted spin-flip excitations",/,2x,35("="))')
case(3)
do ist = 1, nstates
do jst = ist, nstates
call get_mrsf_transition_density(infos, trden(:,:,ist,jst), bvec_mo, ist, jst)
end do
end do
squared_S(:) = 2.0_dp
call get_mrsf_transitions(trans, nocca, noccb, nbf)
write(*,'(/,2x,35("="),/,2x,&
&"Spin-adapted spin-flip excitations",/,2x,35("="))')
case(5)
call get_transition_density(trden, bvec_mo, nbf, noccb, nocca, nstates)
squared_S(:) = 6.0_dp
call get_transitions(trans, noccb, nocca, nbf)
write(*,'(/,2x,35("="),/,2x,&
&"Beta -> Alpha spin-flip excitations",/,2x,35("="))')
case default
error stop "Unknown mrst value"
end select
call get_transition_dipole(basis, dip, mo_a, trden, nstates)
mrsf_energies = eex(1:nstates)
bvec_mo_out = bvec_mo(:,1:nstates)
infos%mol_energy%excited_energy = mrsf_energies(infos%tddft%target_state)
call print_results(infos, bvec_mo, eex, trans, dip, squared_S, nstates)
call flush(iw)
call int2_driver%clean()
call measure_time(print_total=1, log_unit=iw)
close(iw)
end subroutine tdhf_mrsf_energy