TABLE OF CONTENTS
- ABINIT/m_skw
- m_skw/find_rstar_gen
- m_skw/mkstar
- m_skw/mkstar_dk1
- m_skw/mkstar_dk2
- m_skw/skw_eval_bks
- m_skw/skw_free
- m_skw/skw_ncwrite
- m_skw/skw_new
- m_skw/skw_print
- m_skw/skw_t
ABINIT/m_skw [ 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 ]
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