genvchi0u.f90

Go to the documentation of this file.

! Copyright (C) 2026 Wenhan Chen, J. K. Dewhurst and S. Sharma.
! This file is distributed under the terms of the GNU General Public License.
! See the file COPYING for license details.

subroutine genvchi0u(ik0,lock,pmat,evecu,vchi0)
use modmain
use modulr
use modomp
implicit none
! arguments
integer, intent(in) :: ik0
integer(omp_lock_kind), intent(inout) :: lock(nwrf)
complex(8), intent(in) :: pmat(nstsv,nstsv,3),evecu(nstulr,nstulr)
complex(4), intent(inout) :: vchi0(nqpt+2,nqpt+2,nwrf)
! local variables
integer nm,nm2,ist,jst,iw,nthd
integer ikpa,jkpa,iq,jq,i,j,k
integer i1,i2,i3,j1,j2,j3
real(8) ei,ej,eij,t1
complex(8) pmatu(3),a(3,3),z1,z2
! automatic arrays
complex(8) zrq(nqpt),y(nstsv)
complex(4) cw(nwrf),b(-1:nqpt,-1:nqpt)
! external functions
complex(8), external :: zdotc
! matrix size
nm=nqpt+2
nm2=nm**2
call holdthd(nstulr,nthd)
!$OMP PARALLEL DO DEFAULT(SHARED) &
!$OMP PRIVATE(zrq,y,cw,a,b,pmatu) &
!$OMP PRIVATE(jst,ei,ej,eij,t1) &
!$OMP PRIVATE(ikpa,jkpa,i1,i2,i3,j1,j2,j3) &
!$OMP PRIVATE(iq,jq,i,j,k,z1,z2,iw) &
!$OMP SCHEDULE(DYNAMIC) &
!$OMP NUM_THREADS(nthd)
do ist=1,nstulr
  ei=evalu(ist,ik0)
  do jst=1,nstulr
    t1=(wkptnr/(nkpa*omega))*(occulr(ist,ik0)-occulr(jst,ik0))
    if (abs(t1) < 1.d-8) cycle
    ej=evalu(jst,ik0)
    eij=ei-ej
! frequency-dependent part in response function formula for all frequencies
    cw(1:nwrf)=t1/(eij+wrf(1:nwrf))
! determine ρ(Q) = ⟨α|exp(-iQ⋅r)|β⟩ and P = ⟨α|-i∇|β⟩
    zrq(1:nqpt)=0.d0
    pmatu(1:3)=0.d0
    do jkpa=1,nkpa
      j1=ivq(1,jkpa); j2=ivq(2,jkpa); j3=ivq(3,jkpa)
      do ikpa=1,nkpa
        i1=j1-ivq(1,ikpa); i2=j2-ivq(2,ikpa); i3=j3-ivq(3,ikpa)
        iq=ivqiq(i1,i2,i3)
        i=(ikpa-1)*nstsv+1
        j=(jkpa-1)*nstsv+1
        z1=zdotc(nstsv,evecu(i,ist),1,evecu(j,jst),1)
        zrq(iq)=zrq(iq)+z1
        if (ikpa == jkpa) then
          do k=1,3
            call zgemv('N',nstsv,nstsv,zone,pmat(:,:,k),nstsv,evecu(i,jst),1, &
             zzero,y,1)
            z2=zdotc(nstsv,evecu(i,ist),1,y,1)
            pmatu(k)=pmatu(k)+z1*vqc(k,ikpa)+z2
          end do
        end if
      end do
    end do
! Hermitian part of body
    do jq=2,nqpt
      b(2:jq-1,jq)=conjg(b(jq,2:jq-1))
      z1=gclq(jq)*conjg(zrq(jq))
      b(jq:nqpt,jq)=gclq(jq:nqpt)*zrq(jq:nqpt)*z1
    end do
    b(-1:1,-1:1)=0.e0
! calculate 3 × ngrf wings of matrix
    t1=-sqrt(fourpi)/eij
    do i=-1,1
      z1=t1*pmatu(i+2)
      b(i,2:nqpt)=z1*conjg(zrq(2:nqpt))*gclq(2:nqpt)
      do j=2,nqpt
        b(j,i)=conjg(b(i,j))
      end do
    end do
! add to body and wings of the response function
    do iw=1,nwrf
      call omp_set_lock(lock(iw))
      call caxpy(nm2,cw(iw),b,1,vchi0(1,1,iw),1)
      call omp_unset_lock(lock(iw))
    end do
! add to the head of the response function
    t1=fourpi/eij**2
    do j=1,3
      z1=t1*conjg(pmatu(j))
      a(1:3,j)=z1*pmatu(1:3)
    end do
    do iw=1,nwrf
      call omp_set_lock(lock(iw))
      vchi0(1:3,1:3,iw)=vchi0(1:3,1:3,iw)+cw(iw)*a(1:3,1:3)
      call omp_unset_lock(lock(iw))
    end do
  end do
end do
!$OMP END PARALLEL DO
call freethd(nthd)
end subroutine