TABLE OF CONTENTS


ABINIT/m_skw [ Modules ]

[ Top ] [ Modules ]

NAME

  m_skw

FUNCTION

  Shankland-Koelling-Wood Fourier interpolation scheme.

COPYRIGHT

 Copyright (C) 2008-2022 ABINIT group (MG)
 This file is distributed under the terms of the
 GNU General Public License, see ~abinit/COPYING
 or http://www.gnu.org/copyleft/gpl.txt .

SOURCE

16 #if defined HAVE_CONFIG_H
17 #include "config.h"
18 #endif
19 
20 #include "abi_common.h"
21 
22 module m_skw
23 
24  use defs_basis
25  use m_errors
26  use m_abicore
27  use m_xmpi
28  use m_crystal
29  use m_sort
30  use m_nctk
31 #ifdef HAVE_NETCDF
32  use netcdf
33 #endif
34 
35  use m_fstrings,       only : itoa, sjoin, ktoa, yesno, ftoa
36  use m_special_funcs,  only : abi_derfc
37  use m_time,           only : cwtime, cwtime_report
38  use m_numeric_tools,  only : imax_loc, vdiff_t, vdiff_eval, vdiff_print
39  use m_bz_mesh,        only : isamek
40  use m_gsphere,        only : get_irredg
41 
42  implicit none
43 
44  private

m_skw/find_rstar_gen [ Functions ]

[ Top ] [ m_skw ] [ Functions ]

NAME

  find_rstar_gen

FUNCTION

  Find the R-space points generating the stars.
  Set skw%nr and skw%rpts.

INPUTS

  cryst<crystal_t>=Crystalline structure.
  nrwant=Number of R-space points wanted
  rmax(3)=Max reduced components of supercell.
  comm=MPI communicator.

OUTPUT

  or2vals(skw%nr)=||R||**2

SOURCE

796 subroutine find_rstar_gen(skw, cryst, nrwant, rmax, or2vals, comm)
797 
798 !Arguments ------------------------------------
799 !scalars
800  type(skw_t),intent(inout) :: skw
801  type(crystal_t),intent(in) :: cryst
802  integer,intent(in) :: nrwant,comm
803 !arrays
804  integer,intent(in) :: rmax(3)
805  real(dp),allocatable,intent(out) :: or2vals(:)
806 
807 !Local variables-------------------------------
808 !scalars
809  integer :: cnt,nstars,i1,i2,i3,msize,ir,nsh,ish,ss,ee,nst,ierr,nprocs,my_rank,ii
810  real(dp) :: r2_prev
811  !character(len=500) :: msg
812 !arrays
813  integer,allocatable :: iperm(:),rtmp(:,:),rgen(:,:),r2sh(:),shlim(:),sh_start(:),sh_stop(:)
814  integer,allocatable :: recvcounts(:),displs(:),recvbuf(:,:)
815  real(dp),allocatable :: r2tmp(:),cnorm(:)
816 
817 ! *********************************************************************
818 
819  nprocs = xmpi_comm_size(comm); my_rank = xmpi_comm_rank(comm)
820 
821  msize = product(2*rmax + 1)
822  ABI_MALLOC(rtmp, (3, msize))
823  ABI_MALLOC(r2tmp, (msize))
824 
825  cnt = 0
826  do i3=-rmax(3),rmax(3)
827    do i2=-rmax(2),rmax(2)
828      do i1=-rmax(1),rmax(1)
829        cnt = cnt + 1
830        rtmp(:, cnt) = [i1,i2,i3]
831        r2tmp(cnt) = dot_product(rtmp(:,cnt), matmul(cryst%rmet, rtmp(:,cnt)))
832      end do
833    end do
834  end do
835 
836  ! Sort r2tmp
837  ABI_MALLOC(iperm, (msize))
838  iperm = [(i1, i1=1,msize)]
839  call sort_dp(msize, r2tmp, iperm, tol12)
840 
841  ! Find R-points generating the stars.
842  ABI_MALLOC(rgen, (3, msize))
843  do ir=1,msize
844    rgen(:,ir) = rtmp(:,iperm(ir))
845  end do
846  rtmp = rgen
847  ABI_FREE(iperm)
848 
849  ABI_MALLOC(r2sh, (msize))     ! Correspondence between R and the shell index.
850  ABI_MALLOC(shlim, (msize+1))  ! For each shell, the index of the initial G-vector.
851  nsh = 1; r2sh(1) = 1; shlim(1) = 1; r2_prev = zero
852  do ir=2,msize
853    if (abs(r2tmp(ir) - r2_prev) > r2tmp(ir) * tol8) then
854      r2_prev = r2tmp(ir); nsh = nsh + 1; shlim(nsh) = ir
855      !write(std_out,*)"nsh: ",shlim(nsh) - shlim(nsh-1)
856    end if
857    r2sh(ir) = nsh
858  end do
859  shlim(nsh+1) = msize + 1
860  ABI_FREE(r2tmp)
861  ABI_FREE(r2sh)
862 
863  !call get_irredg(msize, skw%ptg_nsym, +1, cryst%rprimd, skw%ptg_symrel, rtmp, nstars, rgen, cnorm)
864  !write(66,*)nstars; do ish=1,nstars; write(66,*)rgen(:,ish); end do
865 
866  ! Distribute shells among processor so that we can parallelize the search algorithm.
867  ! Each proc works on a contigous block of shells, then we have to gather the results.
868  ABI_MALLOC(sh_start, (0:nprocs-1))
869  ABI_MALLOC(sh_stop, (0:nprocs-1))
870  call xmpi_split_work2_i4b(nsh, nprocs, sh_start, sh_stop)
871 
872  ABI_MALLOC(cnorm, (msize))
873  nstars = 0
874  do ish=sh_start(my_rank),sh_stop(my_rank)
875    ss = shlim(ish); ee = shlim(ish+1) - 1; msize = ee - ss + 1
876    call get_irredg(msize, skw%ptg_nsym, + 1, cryst%rprimd, skw%ptg_symrel, rtmp(:,ss:), &
877      nst, rgen(:,nstars+1:), cnorm(nstars+1:))
878    nstars = nstars + nst
879  end do
880 
881  ABI_FREE(cnorm)
882  ABI_FREE(sh_start)
883  ABI_FREE(sh_stop)
884  ABI_FREE(rtmp)
885  ABI_FREE(shlim)
886 
887  if (nprocs > 1) then
888    ! Collect star functions.
889    ABI_MALLOC(recvcounts, (nprocs))
890    recvcounts = 0; recvcounts(my_rank+1) = 3 * nstars
891    call xmpi_sum(recvcounts, comm, ierr)
892    ABI_MALLOC(displs, (nprocs))
893    displs(1) = 0
894    do ii=2,nprocs
895      displs(ii) = sum(recvcounts(:ii-1))
896    end do
897    call xmpi_sum(nstars, nst, comm, ierr)   ! Now nst is the total number of star functions.
898    ABI_MALLOC(recvbuf, (3, nst))
899    call xmpi_allgatherv(rgen, 3*nstars, recvbuf, recvcounts, displs, comm, ierr)
900    ABI_FREE(recvcounts)
901    ABI_FREE(displs)
902    nstars = nst
903    rgen(:,1:nstars) = recvbuf
904    ABI_FREE(recvbuf)
905  end if
906  !if (my_rank == 0) then
907  !  write(67,*)"nstars",nstars,"nsh",nsh; do ish=1,nstars; write(67,*)rgen(:,ish); end do
908  !end if
909 
910  ! Store rpts and compute ||R||**2.
911  skw%nr = min(nstars, nrwant)
912  if (allocated(skw%rpts)) then
913    ABI_FREE(skw%rpts)
914  end if
915  ABI_MALLOC(skw%rpts, (3, skw%nr))
916  skw%rpts = rgen(:,1:skw%nr)
917  ABI_MALLOC(or2vals, (skw%nr))
918  do ir=1,skw%nr
919    or2vals(ir) = dot_product(skw%rpts(:,ir), matmul(cryst%rmet, skw%rpts(:,ir)))
920  end do
921 
922  ABI_FREE(rgen)
923 
924 end subroutine find_rstar_gen

m_skw/mkstar [ Functions ]

[ Top ] [ m_skw ] [ Functions ]

NAME

  mkstar

FUNCTION

  Compute the star function for k-point kpt

INPUTS

  kpt(3)=K-point in reduced coordinates.

OUTPUT

  srk(%nr)=Star function for this k-point.

SOURCE

635 subroutine mkstar(skw, kpt, srk)
636 
637 !Arguments ------------------------------------
638 !scalars
639  type(skw_t),intent(in) :: skw
640 !arrays
641  real(dp),intent(in) :: kpt(3)
642  complex(dpc),intent(out) :: srk(skw%nr)
643 
644 !Local variables-------------------------------
645 !scalars
646  integer :: ir,isym
647 !arrays
648  real(dp) :: sk(3)
649 
650 ! *********************************************************************
651 
652  srk = zero
653  do isym=1,skw%ptg_nsym
654    sk = two_pi * matmul(transpose(skw%ptg_symrel(:,:,isym)), kpt)
655    do ir=1,skw%nr
656      srk(ir) = srk(ir) + exp(j_dpc * dot_product(sk, skw%rpts(:,ir)))
657    end do
658  end do
659  srk = srk / skw%ptg_nsym
660 
661 end subroutine mkstar

m_skw/mkstar_dk1 [ Functions ]

[ Top ] [ m_skw ] [ Functions ]

NAME

  mkstar_dk1

FUNCTION

  Compute the 1st derivative of the star function wrt k

INPUTS

  kpt(3)=K-point in reduced coordinates.

OUTPUT

  srk_dk1(%nr,3)=Derivative of the star function wrt k in reduced coordinates.

SOURCE

681 subroutine mkstar_dk1(skw, kpt, srk_dk1)
682 
683 !Arguments ------------------------------------
684 !scalars
685  type(skw_t),intent(in) :: skw
686 !arrays
687  real(dp),intent(in) :: kpt(3)
688  complex(dpc),intent(out) :: srk_dk1(skw%nr,3)
689 
690 !Local variables-------------------------------
691 !scalars
692  integer :: ir,isym
693 !arrays
694  real(dp) :: sk(3)
695  complex(dpc) :: work(3,skw%nr)
696 
697 ! *********************************************************************
698 
699  work = zero
700  do isym=1,skw%ptg_nsym
701    sk = two_pi * matmul(transpose(skw%ptg_symrel(:,:,isym)), kpt)
702    do ir=1,skw%nr
703      work(:,ir) = work(:,ir) + exp(j_dpc * dot_product(sk, skw%rpts(:,ir))) * &
704         matmul(skw%ptg_symrel(:,:,isym), skw%rpts(:,ir))
705    end do
706  end do
707  work = j_dpc * work / skw%ptg_nsym
708  srk_dk1 = transpose(work)
709 
710 end subroutine mkstar_dk1

m_skw/mkstar_dk2 [ Functions ]

[ Top ] [ m_skw ] [ Functions ]

NAME

  mkstar_dk2

FUNCTION

  Compute the 2st derivatives of the star function wrt k

INPUTS

  kpt(3)=K-point in reduced coordinates.

OUTPUT

  srk_dk2(%nr,3,3)=2nd derivatives of the star function wrt k in reduced coordinates.

SOURCE

730 subroutine mkstar_dk2(skw, kpt, srk_dk2)
731 
732 !Arguments ------------------------------------
733 !scalars
734  type(skw_t),intent(in) :: skw
735 !arrays
736  real(dp),intent(in) :: kpt(3)
737  complex(dpc),intent(out) :: srk_dk2(skw%nr,3,3)
738 
739 !Local variables-------------------------------
740 !scalars
741  integer :: ir,isym,ii,jj
742  complex(dpc) :: eiskr
743 !arrays
744  integer :: sr(3)
745  real(dp) :: sk(3)
746  complex(dpc) :: work(3,3,skw%nr)
747 
748 ! *********************************************************************
749 
750  work = zero
751  do isym=1,skw%ptg_nsym
752    sk = two_pi * matmul(transpose(skw%ptg_symrel(:,:,isym)), kpt)
753    do ir=1,skw%nr
754      sr = matmul(skw%ptg_symrel(:,:,isym), skw%rpts(:,ir))
755      eiskr = exp(j_dpc * dot_product(sk, skw%rpts(:,ir)))
756      do jj=1,3
757        do ii=1,jj
758          work(ii,jj,ir) = work(ii,jj,ir) + eiskr * sr(ii) * sr(jj)
759        end do
760      end do
761    end do
762  end do
763  work = - work / skw%ptg_nsym
764 
765  do jj=1,3
766    do ii=1,jj
767      srk_dk2(:, ii, jj) = work(ii, jj, :)
768      if (ii /= jj) srk_dk2(:,jj,ii) = work(:,ii,jj)
769    end do
770  end do
771 
772 end subroutine mkstar_dk2

m_skw/skw_eval_bks [ Functions ]

[ Top ] [ m_skw ] [ Functions ]

NAME

  skw_eval_bks

FUNCTION

  Interpolate the energies for an arbitrary k-point and spin with slow FT.

INPUTS

  band=Band index (global index associated to the input eigenvalues, i.e. independent of band_block)
  kpt(3)=K-point in reduced coordinates.
  spin=Spin index.

OUTPUT

  oeig=interpolated eigenvalues
    Note that oeig is not necessarily sorted in ascending order.
    The routine does not reorder the interpolated eigenvalues
    to be consistent with the interpolation of the derivatives.
  [oder1(3)]=First-order derivatives wrt k in reduced coordinates.
  [oder2(3,3)]=Second-order derivatives wrt k in reduced coordinates.

SOURCE

525 subroutine skw_eval_bks(skw, band, kpt, spin, oeig, oder1, oder2)
526 
527 !Arguments ------------------------------------
528 !scalars
529  integer,intent(in) :: band,spin
530  class(skw_t),intent(inout) :: skw
531 !arrays
532  real(dp),intent(in) :: kpt(3)
533  real(dp),intent(out) :: oeig
534  real(dp),optional,intent(out) :: oder1(3),oder2(3,3)
535 
536 !Local variables-------------------------------
537 !scalars
538  integer :: ii,jj,ib
539 ! *********************************************************************
540 
541  ib = band - skw%band_block(1) + 1
542  ABI_CHECK(ib >= 1 .and. ib <= skw%bcount, sjoin("out of range band:", itoa(band)))
543 
544  ! Compute star function for this k-point (if not already in memory)
545  if (any(kpt /= skw%cached_kpt)) then
546    call mkstar(skw, kpt, skw%cached_srk)
547    skw%cached_kpt = kpt
548  end if
549 
550  oeig = dot_product(conjg(skw%coefs(:,ib,spin)), skw%cached_srk)
551 
552  ! TODO: Test Derivatives
553  if (present(oder1)) then
554    ! Compute first-order derivatives.
555    if (any(kpt /= skw%cached_kpt_dk1)) then
556      call mkstar_dk1(skw, kpt, skw%cached_srk_dk1)
557      skw%cached_kpt_dk1 = kpt
558    end if
559 
560    do ii=1,3
561      oder1(ii) = dot_product(conjg(skw%coefs(:,ib,spin)), skw%cached_srk_dk1(:,ii)) * two_pi
562    end do
563  end if
564 
565  if (present(oder2)) then
566    ! Compute second-order derivatives.
567    if (any(kpt /= skw%cached_kpt_dk2)) then
568      call mkstar_dk2(skw, kpt, skw%cached_srk_dk2)
569      skw%cached_kpt_dk2 = kpt
570    end if
571 
572    oder2 = zero
573    do jj=1,3
574      do ii=1,jj
575        oder2(ii, jj) = dot_product(conjg(skw%coefs(:,ib,spin)), skw%cached_srk_dk2(:,ii,jj)) * two_pi**2
576        if (ii /= jj) oder2(jj, ii) = oder2(ii, jj)
577      end do
578    end do
579  end if
580 
581 end subroutine skw_eval_bks

m_skw/skw_free [ Functions ]

[ Top ] [ m_skw ] [ Functions ]

NAME

  skw_free

FUNCTION

  Free memory

SOURCE

595 subroutine skw_free(skw)
596 
597 !Arguments ------------------------------------
598 !scalars
599  class(skw_t),intent(inout) :: skw
600 
601 ! *********************************************************************
602 
603  ABI_SFREE(skw%rpts)
604  ABI_SFREE(skw%ptg_symrel)
605  ABI_SFREE(skw%ptg_symrec)
606  ABI_SFREE(skw%coefs)
607 
608  ABI_SFREE(skw%cached_srk)
609  skw%cached_kpt = huge(one)
610  ABI_SFREE(skw%cached_srk_dk1)
611  skw%cached_kpt_dk1 = huge(one)
612  ABI_SFREE(skw%cached_srk_dk2)
613  skw%cached_kpt_dk2 = huge(one)
614 
615 end subroutine skw_free

m_skw/skw_ncwrite [ Functions ]

[ Top ] [ m_skw ] [ Functions ]

NAME

 skw_ncwrite

FUNCTION

   Write the object in netcdf format

INPUTS

  ncid=NC file handle.
  [prefix]=String prepended to netcdf dimensions/variables (HDF5 poor-man groups)
   "skw" if not specified.

OUTPUT

  Only writing

SOURCE

448 integer function skw_ncwrite(self, ncid, prefix) result(ncerr)
449 
450 !Arguments ------------------------------------
451 !scalars
452  class(skw_t),intent(in) :: self
453  integer,intent(in) :: ncid
454  character(len=*),optional,intent(in) :: prefix
455 
456 #ifdef HAVE_NETCDF
457 !Local variables-------------------------------
458 !scalars
459  character(len=500) :: prefix_
460 !arrays
461  real(dp),allocatable :: real_coefs(:,:,:,:)
462 
463 ! *************************************************************************
464 
465  prefix_ = "skw"; if (present(prefix)) prefix_ = trim(prefix)
466 
467  ! Define dimensions.
468  ncerr = nctk_def_dims(ncid, [ &
469    nctkdim_t("nr", self%nr), nctkdim_t("nkpt", self%nkpt), nctkdim_t("bcount", self%bcount), &
470    nctkdim_t("nsppol", self%nsppol)], &
471    defmode=.True., prefix=prefix_)
472  NCF_CHECK(ncerr)
473 
474  ncerr = nctk_def_arrays(ncid, [ &
475   ! Atomic structure and symmetry operations
476   nctkarr_t("rpts", "dp", "three, number_of_cartesian_directions, number_of_vectors"), &
477   nctkarr_t("kpts", "dp", "three, nkpt"), &
478   nctkarr_t("coefs", "dp", "two, nr, bcount, nsppol") &
479  ], prefix=prefix_)
480  NCF_CHECK(ncerr)
481 
482  ! Write data.
483  NCF_CHECK(nctk_set_datamode(ncid))
484  NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, pre("rpts")), self%rpts))
485  !NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, pre("kpts")), self%kpts))
486  ABI_MALLOC(real_coefs, (2, self%nr, self%bcount, self%nsppol))
487  real_coefs(1,:,:,:) = real(self%coefs); real_coefs(2,:,:,:) = aimag(self%coefs)
488  NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, pre("coefs")), real_coefs))
489  ABI_FREE(real_coefs)
490 
491 contains
492   pure function pre(istr) result(ostr)
493     character(len=*),intent(in) :: istr
494     character(len=len_trim(prefix_) + len_trim(istr)+1) :: ostr
495     ostr = trim(prefix_) // trim(istr)
496   end function pre
497 
498 #endif
499 
500 end function skw_ncwrite

m_skw/skw_new [ Functions ]

[ Top ] [ m_skw ] [ Functions ]

NAME

  skw_new

FUNCTION

  Initialize the object.

INPUTS

  cryst<crystal_t>=Crystalline structure.
  params(:)
     params(1): Ratio between star functions and ab-initio k-points.
     params(2:3): Activate Fourier filtering (Eq 9 of PhysRevB.61.1639 [[cite:Uehara2000]]) if params(2) > tol6
       params(2)=rcut, params(3) = rsigma
  cplex=1 if time reversal can be used, 2 otherwise.
  nband=Total Number of bands in the eig array.
  nkpt=Number of ab-initio k-points.
  nsppol=Number of independent spin polarizations.
  kpts(3,nkpt)=ab-initio k-points in reduced coordinates.
  eig(nband,nkpt,nsppol)=ab-initio eigenvalues.
  band_block(2)=Initial and final band index to interpolate. If [0,0], all bands are used
    This is a global variable i.e. all MPI procs MUST call the routine with the same value.
  comm=MPI communicator

SOURCE

166 type(skw_t) function skw_new(cryst, params, cplex, nband, nkpt, nsppol, kpts, eig, band_block, comm) result(new)
167 
168 !Arguments ------------------------------------
169 !scalars
170  integer,intent(in) :: cplex,nband,nkpt,nsppol,comm
171  real(dp),intent(in) :: params(:)
172  type(crystal_t),intent(in) :: cryst
173 !arrays
174  integer,intent(in) :: band_block(2)
175  real(dp),intent(in) :: kpts(3,nkpt)
176  real(dp),intent(in) :: eig(nband,nkpt,nsppol)
177 
178 !Local variables-------------------------------
179 !scalars
180  integer,parameter :: master=0,prtvol=1
181  integer :: my_rank,nprocs,cnt,bstop,bstart,bcount,lwork
182  integer :: ir,ik,ib,ii,jj,nr,band,spin,ierr,lpratio,nrwant
183  real(dp),parameter :: c1=0.25_dp,c2=0.25_dp
184  real(dp) :: r2,r2min,mare,mae_meV,adiff_meV,rel_err,rcut,rsigma
185  real(dp) :: cpu_tot,wall_tot,gflops_tot,cpu,wall,gflops,rval
186  character(len=500) :: fmt,msg
187 !arrays
188  integer :: rmax(3)
189  integer,allocatable :: ipiv(:)
190  real(dp) :: list2(2)
191  real(dp),allocatable :: r2vals(:),inv_rhor(:),oeig(:)
192  complex(dpc),allocatable :: srk(:,:),hmat(:,:),lambda(:,:,:),work(:)
193 
194 ! *********************************************************************
195 
196  ABI_CHECK(nkpt > 1, sjoin("nkpt must be > 1 but got:", itoa(nkpt)))
197 
198  call cwtime(cpu_tot, wall_tot, gflops_tot, "start")
199 
200  nprocs = xmpi_comm_size(comm); my_rank = xmpi_comm_rank(comm)
201 
202  ! Get slice of bands to be treated.
203  new%band_block = band_block; if (all(band_block == 0)) new%band_block = [1, nband]
204  bstart = new%band_block(1); bstop = new%band_block(2); bcount = bstop - bstart + 1
205  new%cplex = cplex; new%nkpt = nkpt; new%nsppol = nsppol; new%bcount = bcount
206 
207  ! Get point group operations.
208  call cryst%get_point_group(new%ptg_nsym, new%ptg_symrel, new%ptg_symrec, new%has_inversion, include_timrev=cplex==1)
209 
210  ! -----------------------
211  ! Find nrwant star points
212  ! -----------------------
213  lpratio = int(abs(params(1)))
214  ABI_CHECK(lpratio > 0, "lpratio must be > 0")
215  rmax = nint((one + (lpratio * new%nkpt * new%ptg_nsym) / two) ** third)
216  if (new%has_inversion) then
217    rmax = nint((one + (lpratio * new%nkpt * new%ptg_nsym / 2) / two) ** third)
218  end if
219  nrwant = lpratio * new%nkpt
220 
221  call cwtime(cpu, wall, gflops, "start")
222  do
223    call find_rstar_gen(new, cryst, nrwant, rmax, r2vals, comm)
224    if (new%nr >= nrwant) then
225      !write(std_out,*)"Entered with rmax", rmax," abs(skw%rpts(last)): ", abs(new%rpts(:,new%nr))
226      exit
227    end if
228    write(std_out,*)"rmax: ", rmax," was not large enough to find ", nrwant," R-star points."
229    rmax = 2 * rmax
230    write(std_out,*)"Will try again with enlarged rmax: ",rmax
231    ABI_FREE(r2vals)
232  end do
233  nr = new%nr
234  call cwtime_report(" find_rstar_gen", cpu, wall, gflops)
235 
236  if (my_rank == master) call new%print(std_out)
237 
238  ! Compute (inverse) roughness function.
239  r2min = r2vals(2)
240  ABI_MALLOC(inv_rhor, (nr))
241  do ir=1,nr
242    r2 = r2vals(ir)
243    inv_rhor(ir) = one / ((one - c1 * r2/r2min)**2 + c2 * (r2 / r2min)**3)
244    ! TODO: Test the two versions.
245    !if (params(1) < zero) inv_rhor(ir) = one / (c1 * r2 + c2 * r2**2)
246  end do
247 
248  ! Construct star functions for the ab-initio k-points.
249  ABI_MALLOC(srk, (nr, nkpt))
250  do ik=1,nkpt
251    call mkstar(new, kpts(:,ik), srk(:,ik))
252  end do
253 
254  ! Build H(k,k') matrix (Hermitian)
255  ABI_CALLOC(hmat, (nkpt-1, nkpt-1))
256  cnt = 0
257  do jj=1,nkpt-1
258    do ii=1,jj
259    !do ii=1,nkpt-1
260      cnt = cnt + 1; if (mod(cnt, nprocs) /= my_rank) cycle ! mpi parallelism.
261      do ir=2,nr
262        hmat(ii, jj) = hmat(ii, jj) + &
263          (srk(ir, ii) - srk(ir, nkpt)) * conjg(srk(ir, jj) - srk(ir, nkpt)) * inv_rhor(ir)
264      end do
265    end do
266  end do
267  call xmpi_sum(hmat, comm, ierr)
268 
269  ABI_MALLOC(lambda, (nkpt-1, bcount, nsppol))
270  do spin=1,nsppol
271    do ib=1,bcount
272      band = ib + bstart - 1
273      lambda(:,ib,spin) = eig(band,1:nkpt-1,spin) - eig(band,nkpt,spin)
274    end do
275  end do
276 
277  ! Solve all bands and spins at once [[cite:Pickett1988]]
278  call wrtout(std_out, " Solving system of linear equations to get lambda coeffients (eq. 10 of PRB 38 2721)...", do_flush=.True.)
279  call cwtime(cpu, wall, gflops, "start")
280  ABI_MALLOC(ipiv, (nkpt-1))
281 
282  if (.False.) then
283    ! General complex.
284    call zgesv(nkpt-1, bcount*nsppol, hmat, nkpt-1, ipiv, lambda, nkpt-1, ierr)
285    ABI_CHECK(ierr == 0, sjoin("ZGESV returned:", itoa(ierr)))
286  else
287    ! Hermitian version
288    do ii=1,nkpt-1
289      hmat(ii, ii) = real(hmat(ii, ii))
290    end do
291    lwork = -1
292    ABI_MALLOC(work, (1))
293    call zhesv("U", nkpt-1, bcount*nsppol, hmat, nkpt-1, ipiv, lambda, nkpt-1, work, lwork, ierr)
294    lwork = nint(real(work(1)))
295    ABI_FREE(work)
296    ABI_MALLOC(work, (lwork))
297    call zhesv("U", nkpt-1, bcount*nsppol, hmat, nkpt-1, ipiv, lambda, nkpt-1, work, lwork, ierr)
298    ABI_CHECK(ierr == 0, sjoin("ZHESV returned:", itoa(ierr)))
299    ABI_FREE(work)
300  end if
301  call cwtime_report(" ZHESV", cpu, wall, gflops)
302 
303  ! Compute coefficients
304  ABI_MALLOC(new%coefs, (nr,bcount,nsppol))
305 
306  do spin=1,nsppol
307    do ib=1,bcount
308      band = ib + bstart - 1
309      do ir=2,nr
310        new%coefs(ir,ib,spin) = inv_rhor(ir) * dot_product(srk(ir,:nkpt-1) - srk(ir,nkpt), lambda(:nkpt-1, ib, spin))
311        !new%coefs(ir,ib,spin) = inv_rhor(ir) * dot_product(lambda(:nkpt-1, ib, spin), conjg(srk(ir,:) - srk(ir,nkpt)))
312        !new%coefs(ir,ib,spin) = inv_rhor(ir) * dot_product(lambda(:nkpt-1, ib, spin), conjg(srk(ir,:) - srk(ir,1)))
313      end do
314      new%coefs(1,ib,spin) = eig(band,nkpt,spin) - dot_product(conjg(new%coefs(2:nr, ib,spin)), srk(2:nr, nkpt))
315    end do
316  end do
317 
318  ! Filter high-frequency.
319  if (params(2) > tol6) then
320    rcut = params(2) * sqrt(r2vals(new%nr))
321    rsigma = params(3); if (rsigma <= zero) rsigma = five
322    call wrtout(std_out," Applying filter (Eq 9 of PhysRevB.61.1639)") ! [[cite:Uehara2000]]
323    do ir=2,nr
324      new%coefs(ir,:,:) = new%coefs(ir,:,:) * half * abi_derfc((sqrt(r2vals(ir)) - rcut) / rsigma)
325    end do
326  end if
327 
328  ! Prepare workspace arrays for star functions.
329  new%cached_kpt = huge(one)
330  ABI_MALLOC(new%cached_srk, (new%nr))
331  new%cached_kpt_dk1 = huge(one)
332  ABI_MALLOC(new%cached_srk_dk1, (new%nr, 3))
333  new%cached_kpt_dk2 = huge(one)
334  ABI_MALLOC(new%cached_srk_dk2, (new%nr, 3, 3))
335 
336  ABI_FREE(r2vals)
337  ABI_FREE(srk)
338  ABI_FREE(inv_rhor)
339  ABI_FREE(hmat)
340  ABI_FREE(lambda)
341  ABI_FREE(ipiv)
342 
343  ! Compare ab-initio data with interpolated results.
344  ABI_MALLOC(oeig, (bcount))
345  fmt = sjoin("(a,", itoa(bcount), "(es12.4))")
346  bstop = bstart + bcount - 1
347  mare = zero; mae_meV = zero; cnt = 0
348  call wrtout(std_out, ch10//" Comparing ab-initio energies with SKW interpolated results...")
349  do spin=1,nsppol
350    do ik=1,nkpt
351      cnt = cnt + 1; if (mod(cnt, nprocs) /= my_rank) cycle ! mpi parallelism.
352 
353      do ib=1,bcount
354        band = ib + new%band_block(1) - 1
355        call new%eval_bks(band, kpts(:,ik), spin, oeig(ib))
356 
357        adiff_meV = abs(eig(band,ik,spin) - oeig(ib)); rel_err = zero
358        if (abs(eig(band,ik,spin)) > tol16) rel_err = adiff_meV / abs(eig(band,ik,spin))
359        rel_err = 100 * rel_err; adiff_meV = adiff_meV * Ha_meV
360        mae_meV = mae_meV + adiff_meV; mare = mare + rel_err
361      end do
362 
363      if (prtvol > 0) then
364        ib = imax_loc(eig(bstart:bstop,ik,spin) - oeig)
365        rval = (eig(bstart+ib-1,ik,spin) - oeig(ib)) * Ha_meV
366        write(std_out,"(a,es12.4,2a)") &
367          " SKW maxerr: ", rval, &
368          " (meV), kpt: ", sjoin(ktoa(kpts(:,ik)), "band:",itoa(bstart+ib-1),", spin: ", itoa(spin))
369        !write(std_out,fmt)"-- ref ", eig(bstart:bstop,ik,spin) * Ha_meV
370        !write(std_out,fmt)"-- int ", oeig * Ha_meV
371        !call vdiff_print(vdiff_eval(1, bcount, eig(bstart:bstop,ik,spin), oeig, one))
372      end if
373    end do
374  end do
375  ABI_FREE(oeig)
376 
377  ! Issue warning if error too large.
378  list2 = [mare, mae_meV]; call xmpi_sum(list2, comm, ierr); mare = list2(1); mae_meV = list2(2)
379  cnt = bcount * nkpt * nsppol; mare = mare / cnt; mae_meV = mae_meV / cnt
380  write(std_out,"(2(a,es12.4),a,/)")" MARE: ",mare, ", MAE: ", mae_meV, " (meV)"
381  if (mae_meV > ten) then
382    write(msg,"(2a,2(a,es12.4),a)") &
383      "Large error in SKW interpolation!",ch10," MARE: ",mare, ", MAE: ", mae_meV, " (meV)"
384    call wrtout(ab_out, msg)
385    ABI_WARNING(msg)
386  end if
387 
388  call cwtime_report(" skw_new", cpu_tot, wall_tot, gflops_tot, end_str=ch10)
389 
390 end function skw_new

m_skw/skw_print [ Functions ]

[ Top ] [ m_skw ] [ Functions ]

NAME

  skw_print

FUNCTION

  Print info on object

INPUTS

  unt=Fortran unit number.

OUTPUT

  only writing

SOURCE

410 subroutine skw_print(skw, unt)
411 
412 !Arguments ------------------------------------
413 !scalars
414  class(skw_t),intent(in) :: skw
415  integer,intent(in) :: unt
416 
417 ! *********************************************************************
418 
419  write(unt,"(a)")" === Shankland-Koelling-Wood Fourier interpolation scheme ==="
420  write(unt,"(a)")sjoin(" nsppol", itoa(skw%nsppol), ", cplex:", itoa(skw%cplex))
421  write(unt,"(a)")sjoin(" Number of ab-initio k-points:", itoa(skw%nkpt))
422  write(unt,"(a)")sjoin(" Number of star functions:", itoa(skw%nr))
423  write(unt,"(a)")sjoin(" Stars/Nk ratio:", ftoa(skw%nr * one / skw%nkpt))
424  write(unt,"(a)")sjoin(" Has spatial inversion:", yesno(skw%has_inversion))
425 
426 end subroutine skw_print

m_skw/skw_t [ Types ]

[ Top ] [ m_skw ] [ Types ]

NAME

 skw_t

FUNCTION

  Object implementing the Shankland-Koelling-Wood Fourier interpolation scheme.
  It can be used to interpolate functions in k-space with the periodicity of the
  reciprocal lattice and satisfying F(k) = F(Sk) for each rotation S
  belonging to the point group of the crystal. For readability reason,
  the names of the variables are chosen assuming we are interpolating electronic eigenvalues
  but the same object can be used to interpolate phonons as well. Just use nsppol=1 and nband = 3 * natom

SOURCE

 63  type,public :: skw_t
 64 
 65   integer :: cplex
 66    ! 1 if time-reversal symmetry can be used, 2 otherwise.
 67 
 68   integer :: nr
 69    ! Number of star functions.
 70 
 71   integer :: nkpt
 72    ! Number of ab-initio k-points.
 73 
 74   integer :: ptg_nsym
 75    ! Number of operations in the point group.
 76 
 77   logical :: has_inversion
 78    ! True if the point group contains spatial inversion.
 79 
 80   integer :: band_block(2)
 81    ! Initial and final band index.
 82 
 83   integer :: bcount
 84    ! Number of bands
 85 
 86   integer :: nsppol
 87    ! Number of independent spin polarizations.
 88 
 89   integer,allocatable :: rpts(:,:)
 90    ! rpts(3, nr)
 91    ! Real-space lattice points (in reduced coordinates) ordered with non-decreasing length.
 92 
 93   integer,allocatable :: ptg_symrel(:,:,:)
 94     ! ptg_symrel(3,3,ptg_nsym)
 95     ! operations of the point group (real space).
 96 
 97   integer,allocatable :: ptg_symrec(:,:,:)
 98     ! ptg_symrec(3,3,ptg_nsym)
 99     ! operations of the point group (reciprocal space).
100 
101   complex(dpc),allocatable :: coefs(:,:,:)
102    ! coefs(nr, bcount, nsppol).
103 
104   complex(dpc),allocatable :: cached_srk(:)
105    ! cached_srk(%nr)
106    ! The star function for cached_kpt (used in skw_eval_bks).
107   real(dp) :: cached_kpt(3)
108 
109   complex(dpc),allocatable :: cached_srk_dk1(:,:)
110    ! cached_srk_dk1(%nr, 3)
111    ! The 1d derivative wrt k of the star function for cached_kpt_dk1 (used in skw_eval_bks).
112   real(dp) :: cached_kpt_dk1(3)
113 
114   complex(dpc),allocatable :: cached_srk_dk2(:,:,:)
115    ! cached_srk_dk2(%nr,3,3)
116    ! The 2d derivatives wrt k of the star function for cached_kpt_dk2 (used in skw_eval_bks).
117   real(dp) :: cached_kpt_dk2(3)
118 
119  contains
120 
121    procedure :: print => skw_print
122    ! Print info about object.
123 
124    procedure :: ncwrite => skw_ncwrite
125    ! Write the object in netcdf format
126 
127    procedure :: eval_bks => skw_eval_bks
128    ! Interpolate eigenvalues, 1st, 2nd derivates wrt k, at an arbitrary k-point.
129 
130    procedure :: free => skw_free
131    ! Free memory.
132 
133  end type skw_t