tddftlru.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 tddftlru
use modmain
use modulr
use modmpi
use modomp
implicit none
! local variables
integer ik0,ik,iw,nm
integer ioc,i,j,n,nthd
character(128) fname
! allocatable arrays
integer(omp_lock_kind), allocatable :: lock(:)
complex(8), allocatable :: evecu(:,:),pmat(:,:,:),vchi0(:,:,:)
complex(4), allocatable :: vchi0_sp(:,:,:)
complex(8), allocatable :: eps0(:,:,:),epsi(:,:,:),epsm(:,:,:)
! no k-point reduction
reducek0=reducek
reducek=0
! initialise global variables
call init0
call init1
call init3
! read the regular Kohn-Sham potential from file
call readstate
! generate the first- and second-variational eigenvectors and eigenvalues for
! the k+κ-point set
call genvsig
call linengy
call genapwlofr
call gensocfr
call genevfsv
! initialise the ultra long-range variables
call initulr
! read in the potential STATE_ULR.OUT
call readstulr
! initialise the external Coulomb potential
call vclqinit
! apply required local operations to the potential and magnetic field
call vblocalu
! determine v¹⸍² (used for the symmetric form of v χ₀)
gclq(2:nqpt)=sqrt(gclq(2:nqpt))
! matrix size
nm=nqpt+2
! initialise the OpenMP locks
allocate(lock(nwrf))
do iw=1,nwrf
  call omp_init_lock(lock(iw))
end do
! compute v¹⸍² χ₀ v¹⸍² in single-precision
allocate(vchi0_sp(nm,nm,nwrf))
vchi0_sp(1:nm,1:nm,1:nwrf)=0.e0
! begin parallel loop over k-points
call holdthd(nkpt0/np_mpi,nthd)
!$OMP PARALLEL DEFAULT(SHARED) &
!$OMP PRIVATE(evecu,pmat,ik) &
!$OMP NUM_THREADS(nthd)
allocate(evecu(nstulr,nstulr),pmat(nstsv,nstsv,3))
!$OMP DO SCHEDULE(DYNAMIC)
do ik0=1,nkpt0
! distribute among MPI processes
  if (mod(ik0-1,np_mpi) /= lp_mpi) cycle
!$OMP CRITICAL(tddftlru_)
  write(*,'("Info(tddftlru): ",I0," of ",I0," k-points")') ik0,nkpt0
!$OMP END CRITICAL(tddftlru_)
! central k-point
  ik=(ik0-1)*nkpa+1
! solve the ultra long-range eigenvalue equation
  call eveqnulr(ik0,evecu)
! get the occupation numbers from file
  call getocculr(ik0)
! generate the momentum matrix elements at the central k-point
  call genpmatk(ik,pmat)
! add to v¹⸍² χ₀ v¹⸍²
  call genvchi0u(ik0,lock,pmat,evecu,vchi0_sp)
end do
!$OMP END DO
deallocate(evecu,pmat)
!$OMP END PARALLEL
call freethd(nthd)
! destroy the OpenMP locks
do iw=1,nwrf
  call omp_destroy_lock(lock(iw))
end do
deallocate(lock)
! add vchi0_sp from each process and redistribute
if (np_mpi > 1) then
  n=nm*nm*nwrf
  call mpi_allreduce(mpi_in_place,vchi0_sp,n,mpi_complex,mpi_sum,mpicom,ierror)
end if
! copy to double-precision array
allocate(vchi0(nm,nm,nwrf))
vchi0(1:nm,1:nm,1:nwrf)=vchi0_sp(1:nm,1:nm,1:nwrf)
deallocate(vchi0_sp)
allocate(eps0(nm,nm,nwrf),epsi(nm,nm,nwrf))
! calculate symmetric ϵ₀ = 1 - v¹⸍² χ₀ v¹⸍²
eps0(1:nm,1:nm,1:nwrf)=-vchi0(1:nm,1:nm,1:nwrf)
do i=1,nm
  eps0(i,i,1:nwrf)=eps0(i,i,1:nwrf)+1.d0
end do
! begin loop over frequencies
do iw=1,nwrf
! invert ϵ₀
  call zminv(nm,eps0(:,:,iw))
! left multiply by v¹⸍² χ₀ v¹⸍²
  call zgemm('N','N',nm,nm,nm,zone,vchi0(:,:,iw),nm,eps0(:,:,iw),nm,zzero, &
   epsi(:,:,iw),nm)
! compute ϵ⁻¹ = 1 + v¹⸍² χ v¹⸍²
  do i=1,nm
    epsi(i,i,iw)=1.d0+epsi(i,i,iw)
  end do
end do
if (mp_mpi) then
! find the macroscopic part of ϵ by inverting the 3 × 3 head only
  allocate(epsm(3,3,nwrf))
  do iw=1,nwrf
    epsm(1:3,1:3,iw)=epsi(1:3,1:3,iw)
    call zminv(3,epsm(:,:,iw))
  end do
! write out the macroscopic components
  do ioc=1,noptcomp
    i=optcomp(1,ioc)
    j=optcomp(2,ioc)
    write(fname,'("EPSM_ULR_",2I1,".OUT")') i,j
    open(50,file=trim(fname),form='FORMATTED')
    do iw=2,nwrf
      write(50,'(2G18.10)') wrf(iw)%re,epsm(i,j,iw)%re
    end do
    write(50,*)
    do iw=2,nwrf
      write(50,'(2G18.10)') wrf(iw)%re,epsm(i,j,iw)%im
    end do
    close(50)
  end do
  write(*,*)
  write(*,'("Info(tddftlru):")')
  write(*,'(" Macroscopic dielectric tensor written to EPSM_ULR_ij.OUT")')
  write(*,'(" for components")')
  do ioc=1,noptcomp
    write(*,'("  i = ",I1,", j = ",I1)') optcomp(1:2,ioc)
  end do
  deallocate(epsm)
end if
deallocate(vchi0,eps0,epsi)
! restore original parameters
reducek=reducek0
end subroutine