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_eval_fft
- m_skw/skw_free
- 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-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 ]
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.