TABLE OF CONTENTS


ABINIT/m_skw [ Modules ]

[ Top ] [ Modules ]

NAME

  m_skw

FUNCTION

  Shankland-Koelling-Wood Fourier interpolation scheme.

COPYRIGHT

 Copyright (C) 2008-2018 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 .

PARENTS

CHILDREN

SOURCE

21 #if defined HAVE_CONFIG_H
22 #include "config.h"
23 #endif
24 
25 #include "abi_common.h"
26 
27 module m_skw
28 
29  use defs_basis
30  use m_errors
31  use m_abicore
32  use m_xmpi
33  use m_crystal
34  use m_sort
35 
36  use m_fstrings,       only : itoa, sjoin, ktoa, yesno, ftoa
37  use m_special_funcs,  only : abi_derfc
38  use m_time,           only : cwtime
39  use m_numeric_tools,  only : imax_loc, vdiff_t, vdiff_eval, vdiff_print
40  use m_bz_mesh,        only : isamek
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

PARENTS

      m_skw

CHILDREN

      get_irredg,sort_dp,xmpi_allgatherv,xmpi_split_work2_i4b,xmpi_sum

SOURCE

 983 subroutine find_rstar_gen(skw, cryst, nrwant, rmax, or2vals, comm)
 984 
 985  use m_gsphere,  only : get_irredg
 986 
 987 !This section has been created automatically by the script Abilint (TD).
 988 !Do not modify the following lines by hand.
 989 #undef ABI_FUNC
 990 #define ABI_FUNC 'find_rstar_gen'
 991 !End of the abilint section
 992 
 993  implicit none
 994 
 995 !Arguments ------------------------------------
 996 !scalars
 997  type(skw_t),intent(inout) :: skw
 998  type(crystal_t),intent(in) :: cryst
 999  integer,intent(in) :: nrwant,comm
1000 !arrays
1001  integer,intent(in) :: rmax(3)
1002  real(dp),allocatable,intent(out) :: or2vals(:)
1003 
1004 !Local variables-------------------------------
1005 !scalars
1006  integer :: cnt,nstars,i1,i2,i3,msize,ir,nsh,ish,ss,ee,nst,ierr,nprocs,my_rank,ii
1007  real(dp) :: r2_prev
1008  character(len=500) :: msg
1009 !arrays
1010  integer,allocatable :: iperm(:),rtmp(:,:),rgen(:,:),r2sh(:),shlim(:),sh_start(:),sh_stop(:)
1011  integer,allocatable :: recvcounts(:),displs(:),recvbuf(:,:)
1012  real(dp),allocatable :: r2tmp(:),cnorm(:)
1013 
1014 ! *********************************************************************
1015 
1016  nprocs = xmpi_comm_size(comm); my_rank = xmpi_comm_rank(comm)
1017 
1018  msize = product(2*rmax + 1)
1019  ABI_MALLOC(rtmp, (3, msize))
1020  ABI_MALLOC(r2tmp, (msize))
1021 
1022  cnt = 0
1023  do i3=-rmax(3),rmax(3)
1024    do i2=-rmax(2),rmax(2)
1025      do i1=-rmax(1),rmax(1)
1026        cnt = cnt + 1
1027        rtmp(:, cnt) = [i1,i2,i3]
1028        r2tmp(cnt) = dot_product(rtmp(:,cnt), matmul(cryst%rmet, rtmp(:,cnt)))
1029      end do
1030    end do
1031  end do
1032 
1033  ! Sort r2tmp
1034  ABI_MALLOC(iperm, (msize))
1035  iperm = [(i1, i1=1,msize)]
1036  call sort_dp(msize, r2tmp, iperm, tol12)
1037 
1038  ! Find R-points generating the stars.
1039  ABI_MALLOC(rgen, (3, msize))
1040  do ir=1,msize
1041    rgen(:,ir) = rtmp(:,iperm(ir))
1042  end do
1043  rtmp = rgen
1044  ABI_FREE(iperm)
1045 
1046  ABI_MALLOC(r2sh, (msize))     ! Correspondence between R and the shell index.
1047  ABI_MALLOC(shlim, (msize+1))  ! For each shell, the index of the initial G-vector.
1048  nsh = 1; r2sh(1) = 1; shlim(1) = 1; r2_prev = zero
1049  do ir=2,msize
1050    if (abs(r2tmp(ir) - r2_prev) > r2tmp(ir) * tol8) then
1051      r2_prev = r2tmp(ir); nsh = nsh + 1; shlim(nsh) = ir
1052      !write(std_out,*)"nsh: ",shlim(nsh) - shlim(nsh-1)
1053    end if
1054    r2sh(ir) = nsh
1055  end do
1056  shlim(nsh+1) = msize + 1
1057  ABI_FREE(r2tmp)
1058  ABI_FREE(r2sh)
1059 
1060  !call get_irredg(msize, skw%ptg_nsym, +1, cryst%rprimd, skw%ptg_symrel, rtmp, nstars, rgen, cnorm)
1061  !write(66,*)nstars; do ish=1,nstars; write(66,*)rgen(:,ish); end do
1062 
1063  ! Distribute shells among processor so that we can parallelize the search algorithm.
1064  ! Each proc works on a contigous block of shells, then we have to gather the results.
1065  ABI_MALLOC(sh_start, (0:nprocs-1))
1066  ABI_MALLOC(sh_stop, (0:nprocs-1))
1067  call xmpi_split_work2_i4b(nsh, nprocs, sh_start, sh_stop, msg, ierr)
1068 
1069  ABI_MALLOC(cnorm, (msize))
1070  nstars = 0
1071  do ish=sh_start(my_rank),sh_stop(my_rank)
1072    ss = shlim(ish); ee = shlim(ish+1) - 1; msize = ee - ss + 1
1073    call get_irredg(msize, skw%ptg_nsym, + 1, cryst%rprimd, skw%ptg_symrel, rtmp(:,ss:), &
1074      nst, rgen(:,nstars+1:), cnorm(nstars+1:))
1075    nstars = nstars + nst
1076  end do
1077 
1078  ABI_FREE(cnorm)
1079  ABI_FREE(sh_start)
1080  ABI_FREE(sh_stop)
1081  ABI_FREE(rtmp)
1082  ABI_FREE(shlim)
1083 
1084  if (nprocs > 1) then
1085    ! Collect star functions.
1086    ABI_MALLOC(recvcounts, (nprocs))
1087    recvcounts = 0; recvcounts(my_rank+1) = 3 * nstars
1088    call xmpi_sum(recvcounts, comm, ierr)
1089    ABI_MALLOC(displs, (nprocs))
1090    displs(1) = 0
1091    do ii=2,nprocs
1092      displs(ii) = sum(recvcounts(:ii-1))
1093    end do
1094    call xmpi_sum(nstars, nst, comm, ierr)   ! Now nst is the total number of star functions.
1095    ABI_MALLOC(recvbuf, (3, nst))
1096    call xmpi_allgatherv(rgen, 3*nstars, recvbuf, recvcounts, displs, comm, ierr)
1097    ABI_FREE(recvcounts)
1098    ABI_FREE(displs)
1099    nstars = nst
1100    rgen(:,1:nstars) = recvbuf
1101    ABI_FREE(recvbuf)
1102  end if
1103  !if (my_rank == 0) then
1104  !  write(67,*)"nstars",nstars,"nsh",nsh; do ish=1,nstars; write(67,*)rgen(:,ish); end do
1105  !end if
1106 
1107  ! Store rpts and compute ||R||**2.
1108  skw%nr = min(nstars, nrwant)
1109  if (allocated(skw%rpts)) then
1110    ABI_FREE(skw%rpts)
1111  end if
1112  ABI_MALLOC(skw%rpts, (3, skw%nr))
1113  skw%rpts = rgen(:,1:skw%nr)
1114  ABI_MALLOC(or2vals, (skw%nr))
1115  do ir=1,skw%nr
1116    or2vals(ir) = dot_product(skw%rpts(:,ir), matmul(cryst%rmet, skw%rpts(:,ir)))
1117  end do
1118 
1119  ABI_FREE(rgen)
1120 
1121 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.

PARENTS

      m_skw

CHILDREN

      get_irredg,sort_dp,xmpi_allgatherv,xmpi_split_work2_i4b,xmpi_sum

SOURCE

777 subroutine mkstar(skw, kpt, srk)
778 
779 
780 !This section has been created automatically by the script Abilint (TD).
781 !Do not modify the following lines by hand.
782 #undef ABI_FUNC
783 #define ABI_FUNC 'mkstar'
784 !End of the abilint section
785 
786  implicit none
787 
788 !Arguments ------------------------------------
789 !scalars
790  type(skw_t),intent(in) :: skw
791 !arrays
792  real(dp),intent(in) :: kpt(3)
793  complex(dpc),intent(out) :: srk(skw%nr)
794 
795 !Local variables-------------------------------
796 !scalars
797  integer :: ir,isym
798 !arrays
799  real(dp) :: sk(3)
800 
801 ! *********************************************************************
802 
803  srk = zero
804  do isym=1,skw%ptg_nsym
805    sk = two_pi * matmul(transpose(skw%ptg_symrel(:,:,isym)), kpt)
806    do ir=1,skw%nr
807      srk(ir) = srk(ir) + exp(j_dpc * dot_product(sk, skw%rpts(:,ir)))
808    end do
809  end do
810  srk = srk / skw%ptg_nsym
811 
812 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.

PARENTS

      m_skw

CHILDREN

      get_irredg,sort_dp,xmpi_allgatherv,xmpi_split_work2_i4b,xmpi_sum

SOURCE

838 subroutine mkstar_dk1(skw, kpt, srk_dk1)
839 
840 
841 !This section has been created automatically by the script Abilint (TD).
842 !Do not modify the following lines by hand.
843 #undef ABI_FUNC
844 #define ABI_FUNC 'mkstar_dk1'
845 !End of the abilint section
846 
847  implicit none
848 
849 !Arguments ------------------------------------
850 !scalars
851  type(skw_t),intent(in) :: skw
852 !arrays
853  real(dp),intent(in) :: kpt(3)
854  complex(dpc),intent(out) :: srk_dk1(skw%nr,3)
855 
856 !Local variables-------------------------------
857 !scalars
858  integer :: ir,isym
859 !arrays
860  real(dp) :: sk(3)
861  complex(dpc) :: work(3,skw%nr)
862 
863 ! *********************************************************************
864 
865  work = zero
866  do isym=1,skw%ptg_nsym
867    sk = two_pi * matmul(transpose(skw%ptg_symrel(:,:,isym)), kpt)
868    do ir=1,skw%nr
869      work(:,ir) = work(:,ir) + exp(j_dpc * dot_product(sk, skw%rpts(:,ir))) * &
870         matmul(skw%ptg_symrel(:,:,isym), skw%rpts(:,ir))
871    end do
872  end do
873  work = j_dpc * work / skw%ptg_nsym
874  srk_dk1 = transpose(work)
875 
876 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.

PARENTS

      m_skw

CHILDREN

      get_irredg,sort_dp,xmpi_allgatherv,xmpi_split_work2_i4b,xmpi_sum

SOURCE

902 subroutine mkstar_dk2(skw, kpt, srk_dk2)
903 
904 
905 !This section has been created automatically by the script Abilint (TD).
906 !Do not modify the following lines by hand.
907 #undef ABI_FUNC
908 #define ABI_FUNC 'mkstar_dk2'
909 !End of the abilint section
910 
911  implicit none
912 
913 !Arguments ------------------------------------
914 !scalars
915  type(skw_t),intent(in) :: skw
916 !arrays
917  real(dp),intent(in) :: kpt(3)
918  complex(dpc),intent(out) :: srk_dk2(skw%nr,3,3)
919 
920 !Local variables-------------------------------
921 !scalars
922  integer :: ir,isym,ii,jj
923  complex(dpc) :: eiskr
924 !arrays
925  integer :: sr(3)
926  real(dp) :: sk(3)
927  complex(dpc) :: work(3,3,skw%nr)
928 
929 ! *********************************************************************
930 
931  work = zero
932  do isym=1,skw%ptg_nsym
933    sk = two_pi * matmul(transpose(skw%ptg_symrel(:,:,isym)), kpt)
934    do ir=1,skw%nr
935      sr = matmul(skw%ptg_symrel(:,:,isym), skw%rpts(:,ir))
936      eiskr = exp(j_dpc * dot_product(sk, skw%rpts(:,ir)))
937      do jj=1,3
938        do ii=1,jj
939          work(ii,jj,ir) = work(ii,jj,ir) + eiskr * sr(ii) * sr(jj)
940        end do
941      end do
942    end do
943  end do
944  work = - work / skw%ptg_nsym
945 
946  do jj=1,3
947    do ii=1,jj
948      srk_dk2(:, ii, jj) = work(ii, jj, :)
949      if (ii /= jj) srk_dk2(:,jj,ii) = work(:,ii,jj)
950    end do
951  end do
952 
953 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.
  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.

PARENTS

      m_ebands,m_ifc,m_skw

CHILDREN

      get_irredg,sort_dp,xmpi_allgatherv,xmpi_split_work2_i4b,xmpi_sum

SOURCE

483 subroutine skw_eval_bks(skw, band, kpt, spin, oeig, oder1, oder2)
484 
485 
486 !This section has been created automatically by the script Abilint (TD).
487 !Do not modify the following lines by hand.
488 #undef ABI_FUNC
489 #define ABI_FUNC 'skw_eval_bks'
490 !End of the abilint section
491 
492  implicit none
493 
494 !Arguments ------------------------------------
495 !scalars
496  integer,intent(in) :: band,spin
497  type(skw_t),intent(inout) :: skw
498 !arrays
499  real(dp),intent(in) :: kpt(3)
500  real(dp),intent(out) :: oeig
501  real(dp),optional,intent(out) :: oder1(3),oder2(3,3)
502 
503 !Local variables-------------------------------
504 !scalars
505  integer :: ii,jj,ib
506 ! *********************************************************************
507 
508  ib = band - skw%band_block(1) + 1
509  !DBG_CHECK(ib >= 1 .and. ib <= skw%bcount, sjoin("out of range band:", itoa(band)))
510 
511  ! Compute star function for this k-point (if not already in memory)
512  if (any(kpt /= skw%cached_kpt)) then
513    call mkstar(skw, kpt, skw%cached_srk)
514    skw%cached_kpt = kpt
515  end if
516 
517  oeig = dot_product(conjg(skw%coefs(:,ib,spin)), skw%cached_srk)
518 
519  ! TODO: Test Derivatives
520  if (present(oder1)) then
521    ! Compute first-order derivatives.
522    if (any(kpt /= skw%cached_kpt_dk1)) then
523      call mkstar_dk1(skw, kpt, skw%cached_srk_dk1)
524      skw%cached_kpt_dk1 = kpt
525    end if
526 
527    do ii=1,3
528      oder1(ii) = dot_product(conjg(skw%coefs(:,ib,spin)), skw%cached_srk_dk1(:,ii))
529    end do
530  end if
531 
532  if (present(oder2)) then
533    ! Compute second-order derivatives.
534    if (any(kpt /= skw%cached_kpt_dk2)) then
535      call mkstar_dk2(skw, kpt, skw%cached_srk_dk2)
536      skw%cached_kpt_dk2 = kpt
537    end if
538 
539    oder2 = zero
540    do jj=1,3
541      do ii=1,jj
542        oder2(ii, jj) = dot_product(conjg(skw%coefs(:,ib,spin)), skw%cached_srk_dk2(:,ii,jj))
543        if (ii /= jj) oder2(jj, ii) = oder2(ii, jj)
544      end do
545    end do
546  end if
547 
548 end subroutine skw_eval_bks

m_skw/skw_eval_fft [ Functions ]

[ Top ] [ m_skw ] [ Functions ]

NAME

  skw_eval_fft

FUNCTION

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

INPUTS

  cryst<crystal_t>=Crystalline structure.
  nfft=Number of points in FFT mesh.
  ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft
  band=Band index.
  spin=Spin index.

OUTPUT

  oeig_mesh(nfft)=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_mesh(3,nfft))]=First-order derivatives wrt k in reduced coordinates.
  [oder2_mesh(3,3,nfft)]=Second-order derivatives wrt k in reduced coordinates.

PARENTS

CHILDREN

      get_irredg,sort_dp,xmpi_allgatherv,xmpi_split_work2_i4b,xmpi_sum

SOURCE

582 subroutine skw_eval_fft(skw, ngfft, nfft, band, spin, oeig_mesh, oder1_mesh, oder2_mesh)
583 
584 
585 !This section has been created automatically by the script Abilint (TD).
586 !Do not modify the following lines by hand.
587 #undef ABI_FUNC
588 #define ABI_FUNC 'skw_eval_fft'
589 !End of the abilint section
590 
591  implicit none
592 
593 !Arguments ------------------------------------
594 !scalars
595  integer,intent(in) :: nfft,band,spin
596  type(skw_t),intent(in) :: skw
597 !arrays
598  integer,intent(in) :: ngfft(18)
599  real(dp),intent(out) :: oeig_mesh(nfft)
600  real(dp),optional,intent(out) :: oder1_mesh(3,nfft)
601  real(dp),optional,intent(out) :: oder2_mesh(3,3,nfft)
602 
603 !Local variables-------------------------------
604 !scalars
605  integer,parameter :: tim_fourdp0=0,paral_kgb0=0
606  integer :: cplex,ix,iy,iz,nx,ny,nz,ldx,ldy,ldz,ifft,ir,ii,jj
607 !arrays
608  real(dp),allocatable :: fofg(:,:),fofr(:),workg(:,:)
609 ! *********************************************************************
610 
611  ! Internal MPI_type needed for calling fourdp!
612  !call initmpi_seq(skw%mpi_enreg)
613  ! Initialize tables to call fourdp in sequential
614  !call init_distribfft_seq(skw%mpi_enreg%distribfft,'c',ngfft(2),ngfft(3),'all')
615  !call init_distribfft_seq(skw%mpi_enreg%distribfft,'f',ngfft(2),ngfft(3),'all')
616 
617  ! Transfer data from the G-sphere to the FFT box.
618  ! Use the following indexing (N means ngfft of the adequate direction)
619  ! 0 1 2 3 ... N/2    -(N-1)/2 ... -1    <= gc
620  ! 1 2 3 4 ....N/2+1  N/2+2    ...  N    <= index ig
621  nx = ngfft(1); ny = ngfft(2); nz = ngfft(3)
622  ldx = ngfft(4); ldy = ngfft(5); ldz = ngfft(6)
623 
624  cplex = 2
625  ABI_MALLOC(fofr, (cplex*nfft))
626  ABI_MALLOC(fofg, (2,nfft))
627 
628  ! TODO: Complete Derivatives
629  ! Decide between fourdp and fftbox API
630  fofg = zero
631  do ir=1,skw%nr
632    ix = skw%rpts(1,ir); if (ix<0) ix=ix+nx; ix=ix+1
633    iy = skw%rpts(2,ir); if (iy<0) iy=iy+ny; iy=iy+1
634    iz = skw%rpts(3,ir); if (iz<0) iz=iz+nz; iz=iz+1
635    ifft = ix + (iy-1)*ldx + (iz-1)*ldx*ldy
636    !band = ib + bstart - 1
637    fofg(1,ifft) = real(skw%coefs(ir, band, spin))
638    fofg(2,ifft) = aimag(skw%coefs(ir, band, spin))
639  end do
640 
641  !call fourdp(cplex, fofg, fofr, +1, skw%mpi_enreg, nfft, ngfft, paral_kgb0, tim_fourdp0)
642  !fofr = fofr / skw%ptg_nsym
643  if (cplex == 1) oeig_mesh = fofr
644  if (cplex == 2) oeig_mesh = fofr(1::2)
645 
646  if (present(oder1_mesh)) then
647    ! Compute first-order derivatives.
648    ABI_MALLOC(workg, (2,nfft))
649    do ii=1,3
650       workg = zero
651       do ir=1,skw%nr
652         !ifft
653         workg(1,ifft) = -fofg(2,ifft) * skw%rpts(ii,ir)
654         workg(2,ifft) =  fofg(1,ifft) * skw%rpts(ii,ir)
655       end do
656       !call fourdp(cplex, workg, fofr, +1, skw%mpi_enreg, nfft, ngfft, paral_kgb0, tim_fourdp0)
657       if (cplex == 1) oder1_mesh(ii,:) = fofr
658       if (cplex == 2) oder1_mesh(ii,:) = fofr(1::2)
659    end do
660    ABI_FREE(workg)
661  end if
662 
663  if (present(oder2_mesh)) then
664    ! Compute second-order derivatives.
665    ABI_MALLOC(workg, (2,nfft))
666    do jj=1,3
667      do ii=1,jj
668        workg = zero
669        do ir=1,skw%nr
670          !ifft
671          workg(:,ifft) = -fofg(:,ifft) * skw%rpts(ii,ir) * skw%rpts(jj,ir)
672        end do
673        !call fourdp(cplex, workg, fofr, +1, skw%mpi_enreg, nfft, ngfft, paral_kgb0, tim_fourdp0)
674        if (cplex == 1) oder2_mesh(ii,jj,:) = fofr
675        if (cplex == 2) oder2_mesh(ii,jj,:) = fofr(1::2)
676        if (ii /= jj) oder2_mesh(jj, ii,:) = oder2_mesh(ii, jj,:)
677      end do
678    end do
679    ABI_FREE(workg)
680  end if
681 
682  ABI_FREE(fofg)
683  ABI_FREE(fofr)
684 
685 end subroutine skw_eval_fft

m_skw/skw_free [ Functions ]

[ Top ] [ m_skw ] [ Functions ]

NAME

  skw_free

FUNCTION

  Free memory

PARENTS

      m_ebands,m_ifc

CHILDREN

      get_irredg,sort_dp,xmpi_allgatherv,xmpi_split_work2_i4b,xmpi_sum

SOURCE

705 subroutine skw_free(skw)
706 
707 
708 !This section has been created automatically by the script Abilint (TD).
709 !Do not modify the following lines by hand.
710 #undef ABI_FUNC
711 #define ABI_FUNC 'skw_free'
712 !End of the abilint section
713 
714  implicit none
715 
716 !Arguments ------------------------------------
717 !scalars
718  type(skw_t),intent(inout) :: skw
719 
720 ! *********************************************************************
721 
722  if (allocated(skw%rpts)) then
723    ABI_FREE(skw%rpts)
724  end if
725  if (allocated(skw%ptg_symrel)) then
726    ABI_FREE(skw%ptg_symrel)
727  end if
728  if (allocated(skw%ptg_symrec)) then
729    ABI_FREE(skw%ptg_symrec)
730  end if
731 
732  if (allocated(skw%coefs)) then
733    ABI_FREE(skw%coefs)
734  end if
735 
736  if (allocated(skw%cached_srk)) then
737    ABI_FREE(skw%cached_srk)
738  end if
739  skw%cached_kpt = huge(one)
740 
741  if (allocated(skw%cached_srk_dk1)) then
742    ABI_FREE(skw%cached_srk_dk1)
743  end if
744  skw%cached_kpt_dk1 = huge(one)
745 
746  if (allocated(skw%cached_srk_dk2)) then
747    ABI_FREE(skw%cached_srk_dk2)
748  end if
749  skw%cached_kpt_dk2 = huge(one)
750 
751 end subroutine skw_free

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

PARENTS

CHILDREN

      sort_dp

SOURCE

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

PARENTS

      m_skw

CHILDREN

      get_irredg,sort_dp,xmpi_allgatherv,xmpi_split_work2_i4b,xmpi_sum

SOURCE

425 subroutine skw_print(skw, unt)
426 
427 
428 !This section has been created automatically by the script Abilint (TD).
429 !Do not modify the following lines by hand.
430 #undef ABI_FUNC
431 #define ABI_FUNC 'skw_print'
432 !End of the abilint section
433 
434  implicit none
435 
436 !Arguments ------------------------------------
437 !scalars
438  type(skw_t),intent(in) :: skw
439  integer,intent(in) :: unt
440 
441 ! *********************************************************************
442 
443  write(unt,"(a)")" === Shankland-Koelling-Wood Fourier interpolation scheme ==="
444  write(unt,"(a)")sjoin(" nsppol", itoa(skw%nsppol), ", cplex:", itoa(skw%cplex))
445  write(unt,"(a)")sjoin(" Number of ab-initio k-points:", itoa(skw%nkpt))
446  write(unt,"(a)")sjoin(" Number of star functions:", itoa(skw%nr))
447  write(unt,"(a)")sjoin(" Stars/Nk ratio:", ftoa(skw%nr * one / skw%nkpt))
448  write(unt,"(a)")sjoin(" Has spatial inversion:", yesno(skw%has_inversion))
449 
450 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 :: nsppol
 84    ! Number of independent spin polarizations.
 85 
 86   integer,allocatable :: rpts(:,:)
 87    ! rpts(3, nr)
 88    ! Real-space lattice points (in reduced coordinates) ordered with non-decreasing length.
 89 
 90   integer,allocatable :: ptg_symrel(:,:,:)
 91     ! ptg_symrel(3,3,ptg_nsym)
 92     ! operations of the point group (real space).
 93 
 94   integer,allocatable :: ptg_symrec(:,:,:)
 95     ! ptg_symrec(3,3,ptg_nsym)
 96     ! operations of the point group (reciprocal space).
 97 
 98   complex(dpc),allocatable :: coefs(:,:,:)
 99    ! coefs(nr, nbcount, nsppol).
100 
101    !type(skcoefs_t),allocatable :: coefs(:,:)
102    ! coefs(mband, 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  end type skw_t
120 
121  public :: skw_new          ! Create new object.
122  public :: skw_print        ! Print info about object.
123  public :: skw_eval_bks     ! Interpolate eigenvalues, 1st, 2nd derivates wrt k, at an arbitrary k-point.
124  public :: skw_free         ! Free memory.