TABLE OF CONTENTS
- ABINIT/m_ddk
- m_ddk/ddk_red2car
- m_ddk/ddkop_apply
- m_ddk/ddkop_free
- m_ddk/ddkop_get_braket
- m_ddk/ddkop_get_vdiag
- m_ddk/ddkop_get_vnondiag
- m_ddk/ddkop_new
- m_ddk/ddkop_setup_spin_kpoint
- m_ddk/ddkop_t
- m_ddk/ddkstore_compute_ddk
- m_ddk/ddkstore_free
- m_ddk/ddkstore_t
- m_ddk/ham_targets_free
ABINIT/m_ddk [ Modules ]
NAME
m_ddk
FUNCTION
Objects and methods to extract data from DDK files. The DDK files are binary (soon also netcdf) files with Hamiltonian derivatives wrt k, and the corresponding wave functions
COPYRIGHT
Copyright (C) 2016-2022 ABINIT group (MJV, HM, 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 . For the initials of contributors, see ~abinit/doc/developers/contributors.txt.
SOURCE
19 #if defined HAVE_CONFIG_H 20 #include "config.h" 21 #endif 22 23 #include "abi_common.h" 24 25 MODULE m_ddk 26 27 use defs_basis 28 use m_abicore 29 use m_errors 30 use m_xmpi 31 use m_nctk 32 use m_hdr 33 use m_dtset 34 use m_krank 35 use m_crystal 36 use m_mpinfo 37 use m_cgtools 38 use m_hamiltonian 39 use m_initylmg 40 use m_ebands 41 use m_pawcprj 42 use m_getgh1c 43 use netcdf 44 45 use m_fstrings, only : strcat, sjoin, itoa, ktoa 46 use m_io_tools, only : iomode_from_fname 47 use m_time, only : cwtime, cwtime_report 48 use defs_abitypes, only : MPI_type 49 use defs_datatypes, only : ebands_t, pseudopotential_type 50 use m_vkbr, only : vkbr_t, nc_ihr_comm, vkbr_init, vkbr_free 51 use m_pawtab, only : pawtab_type 52 use m_wfk, only : wfk_read_ebands !, wfk_read_h1mat 53 use m_wfd, only : wfd_t, wfd_init, wave_t 54 55 implicit none 56 57 private
m_ddk/ddk_red2car [ Functions ]
[ Top ] [ m_ddk ] [ Functions ]
NAME
ddk_red2car
FUNCTION
Convert ddk matrix element from reduced coordinates to cartesian coordinates.
INPUTS
OUTPUT
SOURCE
655 pure subroutine ddk_red2car(rprimd, vred, vcar) 656 657 !Arguments ------------------------------------- 658 real(dp),intent(in) :: rprimd(3,3) 659 real(dp),intent(in) :: vred(2,3) 660 real(dp),intent(out) :: vcar(2,3) 661 662 !Local variables ------------------------------- 663 real(dp) :: vtmp(2,3) 664 665 !************************************************************************ 666 667 ! Go to Cartesian coordinates (same as pmat2cart routine) 668 ! V_cart = 1/(2pi) * Rprimd x V_red 669 ! where V_red is the derivative computed in the DFPT routines (derivative wrt reduced component). 670 vtmp(1,:) = rprimd(:,1)*vred(1,1) & 671 +rprimd(:,2)*vred(1,2) & 672 +rprimd(:,3)*vred(1,3) 673 vtmp(2,:) = rprimd(:,1)*vred(2,1) & 674 +rprimd(:,2)*vred(2,2) & 675 +rprimd(:,3)*vred(2,3) 676 vcar = vtmp / two_pi 677 678 end subroutine ddk_red2car
m_ddk/ddkop_apply [ Functions ]
[ Top ] [ m_ddk ] [ Functions ]
NAME
ddkop_apply
FUNCTION
Apply velocity operator dH/dk to wavefunction in G-space. Store results in object.
INPUTS
eig0nk: Eigenvalue associated to the wavefunction. npw_k: Number of planewaves. nspinor: Number of spinor components. cwave(2,npw_k*nspinor)=input wavefunction in reciprocal space cwaveprj(natom,nspinor*usecprj)=<p_lmn|C> coefficients for wavefunction |C> (and 1st derivatives) if not allocated or size=0, they are locally computed (and not sorted)!!
SIDE EFFECTS
Stores: gh1c(2,npw1*nspinor)= <G|H^(1)|C> or <G|H^(1)-lambda.S^(1)|C> on the k+q sphere (only kinetic+non-local parts if optlocal=0)
SOURCE
875 subroutine ddkop_apply(self, eig0nk, npw_k, nspinor, cwave, cwaveprj) 876 877 !Arguments ------------------------------------ 878 !scalars 879 class(ddkop_t),intent(inout) :: self 880 integer,intent(in) :: npw_k, nspinor 881 real(dp),intent(in) :: eig0nk 882 !arrays 883 real(dp),intent(inout) :: cwave(2,npw_k*nspinor) 884 type(pawcprj_type),intent(inout) :: cwaveprj(:,:) 885 886 !Local variables------------------------------- 887 !scalars 888 integer,parameter :: berryopt0 = 0, optlocal0 = 0, tim_getgh1c = 1, usevnl0 = 0, opt_gvnlx1 = 0 889 integer :: idir, sij_opt, ispinor, ipws, ipw, optnl 890 real(dp) :: eshift 891 !arrays 892 real(dp) :: grad_berry(2,(berryopt0/4)), gvnlx1(2,usevnl0) 893 real(dp),pointer :: dkinpw(:),kinpw1(:) 894 895 !************************************************************************ 896 897 self%eig0nk = eig0nk 898 899 if (self%inclvkb /= 0) then 900 ! optlocal0 = 0: local part of H^(1) is not computed in gh1c=<G|H^(1)|C> 901 ! optnl = 2: non-local part of H^(1) is totally computed in gh1c=<G|H^(1)|C> 902 ! opt_gvnlx1 = option controlling the use of gvnlx1 array: 903 optnl = 2 !; if (self%inclvkb == 0) optnl = 0 904 905 eshift = self%eig0nk - self%dfpt_sciss 906 do idir=1,3 907 sij_opt = self%gs_hamkq(idir)%usepaw 908 call getgh1c(berryopt0, cwave, cwaveprj, self%gh1c(:,:,idir), & 909 grad_berry, self%gs1c(:,:,idir), self%gs_hamkq(idir), gvnlx1, idir, self%ipert, eshift, self%mpi_enreg, optlocal0, & 910 optnl, opt_gvnlx1, self%rf_hamkq(idir), sij_opt, tim_getgh1c, usevnl0) 911 end do 912 913 else 914 ! optnl 0 with DDK does not work as expected. 915 ! So I treat the kinetic term explicitly without calling getgh1c. 916 do idir=1,3 917 kinpw1 => self%gs_hamkq(idir)%kinpw_kp 918 dkinpw => self%rf_hamkq(idir)%dkinpw_k 919 do ispinor=1,nspinor 920 do ipw=1,npw_k 921 ipws = ipw + npw_k*(ispinor-1) 922 if (kinpw1(ipw) < huge(zero)*1.d-11) then 923 self%gh1c(:,ipws,idir) = dkinpw(ipw) * cwave(:,ipws) 924 else 925 self%gh1c(:,ipws,idir) = zero 926 end if 927 end do 928 end do 929 end do 930 end if 931 932 end subroutine ddkop_apply
m_ddk/ddkop_free [ Functions ]
[ Top ] [ m_ddk ] [ Functions ]
NAME
ddkop_free
FUNCTION
Free memory
INPUTS
SOURCE
1103 subroutine ddkop_free(self) 1104 1105 !Arguments ------------------------------------ 1106 !scalars 1107 class(ddkop_t),intent(inout) :: self 1108 1109 !Local variables------------------------------- 1110 !scalars 1111 integer :: idir 1112 1113 !************************************************************************ 1114 1115 ABI_SFREE(self%gh1c) 1116 ABI_SFREE(self%gs1c) 1117 1118 do idir=1,3 1119 call self%gs_hamkq(idir)%free() 1120 call self%htg(idir)%free() 1121 call self%rf_hamkq(idir)%free() 1122 end do 1123 1124 self%mpi_enreg => null() 1125 1126 end subroutine ddkop_free
m_ddk/ddkop_get_braket [ Functions ]
[ Top ] [ m_ddk ] [ Functions ]
NAME
ddkop_get_braket
FUNCTION
Compute diagonal matrix element in Cartesian coordinates.
INPUTS
eig0mk: Eigenvalue associated to the "bra" wavefunction istwkf_k: defines storage of wavefunctions for this k-point npw_k: Number of planewaves. nspinor: Number of spinor components. brag(2,npw_k*nspinor)=input wavefunction in reciprocal space
SOURCE
953 function ddkop_get_braket(self, eig0mk, istwf_k, npw_k, nspinor, brag, mode) result(vk) 954 955 !Arguments ------------------------------------ 956 !scalars 957 class(ddkop_t),intent(in) :: self 958 integer,intent(in) :: istwf_k, npw_k, nspinor 959 real(dp),intent(in) :: eig0mk 960 character(len=*),optional,intent(in) :: mode 961 !arrays 962 real(dp),intent(in) :: brag(2*npw_k*nspinor) 963 real(dp) :: vk(2,3) 964 965 !Local variables------------------------------- 966 !scalars 967 integer :: idir 968 real(dp) :: doti 969 !arrays 970 real(dp) :: dotarr(2), vk_red(2, 3) 971 character(len=50) :: my_mode 972 973 !************************************************************************ 974 975 if (self%usepaw == 0) then 976 ! <u_(iband,k+q)^(0)|H_(k+q,k)^(1)|u_(jband,k)^(0)> (NC psps) 977 do idir=1,3 978 dotarr = cg_zdotc(npw_k * nspinor, brag, self%gh1c(:,:,idir)) 979 if (istwf_k > 1) then 980 !dum = two * j_dpc * AIMAG(dum); if (vkbr%istwfk==2) dum = dum - j_dpc * AIMAG(gamma_term) 981 doti = two * dotarr(2) 982 if (istwf_k == 2 .and. self%mpi_enreg%me_g0 == 1) then 983 ! nspinor always 1 984 ! TODO: Recheck this part but it should be ok. 985 doti = doti - (brag(1) * self%gh1c(2,1,idir) - brag(2) * self%gh1c(1,1,idir)) 986 end if 987 dotarr(2) = doti; dotarr(1) = zero 988 end if 989 vk(:, idir) = dotarr 990 end do 991 else 992 ABI_ERROR("PAW Not Implemented") 993 ! <u_(iband,k+q)^(0)|H_(k+q,k)^(1)-(eig0_k+eig0_k+q)/2.S^(1)|u_(jband,k)^(0)> (PAW) 994 ! eshiftkq = half * (eig0mk - self%eig0nk) 995 ABI_UNUSED(eig0mk) 996 end if 997 998 my_mode = "cart"; if (present(mode)) my_mode = mode 999 select case (mode) 1000 case ("cart") 1001 vk_red = vk 1002 call ddk_red2car(self%rprimd, vk_red, vk) 1003 case ("reduced") 1004 continue 1005 case default 1006 ABI_ERROR(sjoin("Invalid vaue for mode:", mode)) 1007 end select 1008 1009 end function ddkop_get_braket
m_ddk/ddkop_get_vdiag [ Functions ]
[ Top ] [ m_ddk ] [ Functions ]
NAME
ddkop_get_vdiag
FUNCTION
Simplified interface to compute the diagonal matrix element of the velocity operator in cartesian coords.
INPUTS
SOURCE
1025 function ddkop_get_vdiag(self, eig0nk, istwf_k, npw_k, nspinor, cwave, cwaveprj, mode) result(vk) 1026 1027 !Arguments ------------------------------------ 1028 !scalars 1029 class(ddkop_t),intent(inout) :: self 1030 integer,intent(in) :: istwf_k, npw_k, nspinor 1031 real(dp),intent(in) :: eig0nk 1032 character(len=*),optional,intent(in) :: mode 1033 !arrays 1034 real(dp),intent(inout) :: cwave(2,npw_k*nspinor) 1035 type(pawcprj_type),intent(inout) :: cwaveprj(:,:) 1036 real(dp) :: vk(3) 1037 1038 !Local variables------------------------------- 1039 character(len=50) :: my_mode 1040 !arrays 1041 real(dp) :: cvk(2, 3) 1042 1043 !************************************************************************ 1044 1045 my_mode = "cart"; if (present(mode)) my_mode = mode 1046 call self%apply(eig0nk, npw_k, nspinor, cwave, cwaveprj) 1047 cvk = self%get_braket(eig0nk, istwf_k, npw_k, nspinor, cwave, mode=my_mode) 1048 vk = cvk(1, :) 1049 1050 end function ddkop_get_vdiag
m_ddk/ddkop_get_vnondiag [ Functions ]
[ Top ] [ m_ddk ] [ Functions ]
NAME
ddkop_get_vdiag
FUNCTION
Simplified interface to compute the off-diagonal matrix elemente of the velocity operator in cartesian coords.
INPUTS
SOURCE
1064 function ddkop_get_vnondiag(self, eig0nk_bra, istwf_k, npw_k, nspinor, cwave_bra, cwave_ket, cwaveprj, mode) result(cvk) 1065 1066 !Arguments ------------------------------------ 1067 !scalars 1068 class(ddkop_t),intent(inout) :: self 1069 integer,intent(in) :: istwf_k, npw_k, nspinor 1070 real(dp),intent(in) :: eig0nk_bra 1071 character(len=*),optional,intent(in) :: mode 1072 !arrays 1073 real(dp),intent(inout) :: cwave_bra(2,npw_k*nspinor),cwave_ket(2,npw_k*nspinor) 1074 type(pawcprj_type),intent(inout) :: cwaveprj(:,:) 1075 real(dp) :: cvk(2,3) 1076 1077 !Local variables------------------------------- 1078 character(len=50) :: my_mode 1079 !arrays 1080 1081 !************************************************************************ 1082 1083 my_mode = "cart"; if (present(mode)) my_mode = mode 1084 call self%apply(eig0nk_bra, npw_k, nspinor, cwave_ket, cwaveprj) 1085 cvk = self%get_braket(eig0nk_bra, istwf_k, npw_k, nspinor, cwave_bra, mode=my_mode) 1086 1087 end function ddkop_get_vnondiag
m_ddk/ddkop_new [ Functions ]
[ Top ] [ m_ddk ] [ Functions ]
NAME
ddkop_new
FUNCTION
Build new object. Use dtset%inclvkb to determine whether non-local part should be included.
INPUTS
dtset<dataset_type>=All input variables for this dataset. cryst<crystal_t>=Crystal structure. pawtab(ntypat*usepaw)<pawtab_type>=Paw tabulated starting data. psps<pseudopotential_type>=Variables related to pseudopotentials. mpi_enreg=information about MPI parallelization mpw=Maximum number of plane-waves over k-points. ngfft(18)=contain all needed information about 3D FFT
OUTPUT
SOURCE
703 type(ddkop_t) function ddkop_new(dtset, cryst, pawtab, psps, mpi_enreg, mpw, ngfft) result(new) 704 705 !Arguments ------------------------------------ 706 !scalars 707 type(dataset_type),intent(in) :: dtset 708 type(crystal_t),intent(in) :: cryst 709 type(pseudopotential_type),intent(in) :: psps 710 type(MPI_type),target,intent(in) :: mpi_enreg 711 integer,intent(in) :: mpw 712 !arrays 713 integer,intent(in) :: ngfft(18) 714 type(pawtab_type),intent(in) :: pawtab(psps%ntypat*psps%usepaw) 715 716 !Local variables------------------------------- 717 !scalars 718 integer,parameter :: cplex1 = 1 719 integer :: nfft, mgfft, idir 720 721 ! ************************************************************************* 722 723 ABI_CHECK(dtset%usepaw == 0, "PAW not tested/implemented!") 724 725 new%inclvkb = dtset%inclvkb 726 new%usepaw = dtset%usepaw 727 new%ipert = cryst%natom + 1 728 new%dfpt_sciss = dtset%dfpt_sciss 729 new%mpw = mpw 730 731 new%rprimd = cryst%rprimd 732 new%mpi_enreg => mpi_enreg 733 734 ! Not used because vlocal1 is not applied. 735 nfft = product(ngfft(1:3)) 736 mgfft = maxval(ngfft(1:3)) 737 738 ABI_MALLOC(new%gh1c, (2, new%mpw*dtset%nspinor, 3)) 739 ABI_MALLOC(new%gs1c, (2, new%mpw*dtset%nspinor, 3)) 740 741 do idir=1,3 742 ! ==== Initialize most of the Hamiltonian (and derivative) ==== 743 ! 1) Allocate all arrays and initialize quantities that do not depend on k and spin. 744 ! 2) Perform the setup needed for the non-local factors: 745 ! * Norm-conserving: Constant kleimann-Bylander energies are copied from psps to gs_hamk. 746 ! * PAW: Initialize the overlap coefficients and allocate the Dij coefficients. 747 call init_hamiltonian(new%gs_hamkq(idir), psps, pawtab, dtset%nspinor, dtset%nsppol, dtset%nspden, cryst%natom,& 748 cryst%typat, cryst%xred, nfft, mgfft, ngfft, cryst%rprimd, dtset%nloalg) 749 !paw_ij=paw_ij,comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab,mpi_spintab=mpi_enreg%my_isppoltab,& 750 !usecprj=usecprj,ph1d=ph1d,nucdipmom=dtset%nucdipmom,use_gpu_cuda=dtset%use_gpu_cuda) 751 752 ! Prepare application of the NL part. 753 call init_rf_hamiltonian(cplex1, new%gs_hamkq(idir), new%ipert, new%rf_hamkq(idir), has_e1kbsc=.true.) 754 !&paw_ij1=paw_ij1,comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab,& 755 !&mpi_spintab=mpi_enreg%my_isppoltab) 756 end do 757 758 end function ddkop_new
m_ddk/ddkop_setup_spin_kpoint [ Functions ]
[ Top ] [ m_ddk ] [ Functions ]
NAME
ddkop_setup_spin_kpoint
FUNCTION
Prepare internal tables that depend on k-point/spin
INPUTS
dtset<dataset_type>=All input variables for this dataset. cryst<crystal_t>=Crystal structure. psps<pseudopotential_type>=Variables related to pseudopotentials. spin: spin index kpoint(3): K-point in reduced coordinates. istwkf_k: defines storage of wavefunctions for this k-point npw_k: Number of planewaves. kg_k(3,npw_k)=reduced planewave coordinates.
SOURCE
782 subroutine ddkop_setup_spin_kpoint(self, dtset, cryst, psps, spin, kpoint, istwf_k, npw_k, kg_k) 783 784 !Arguments ------------------------------------ 785 !scalars 786 integer,intent(in) :: spin, npw_k, istwf_k 787 class(ddkop_t),intent(inout) :: self 788 type(crystal_t) :: cryst 789 type(dataset_type),intent(in) :: dtset 790 type(pseudopotential_type),intent(in) :: psps 791 !arrays 792 integer,intent(in) :: kg_k(3,npw_k) 793 real(dp),intent(in) :: kpoint(3) 794 795 !Local variables------------------------------- 796 !scalars 797 integer,parameter :: nkpt1=1, nsppol1=1 798 type(mpi_type) :: mpienreg_seq 799 !arrays 800 integer :: npwarr(nkpt1), dummy_nband(nkpt1*nsppol1) 801 integer :: idir, nkpg, nkpg1, useylmgr1, optder !, nylmgr1 802 real(dp),allocatable :: ylm_k(:,:),ylmgr1_k(:,:,:) 803 804 !************************************************************************ 805 806 ABI_CHECK(npw_k <= self%mpw, "npw_k > mpw!") 807 self%kpoint = kpoint 808 809 ! Set up the spherical harmonics (Ylm) at k+q if useylm = 1 810 useylmgr1 = 0; optder = 0 811 if (psps%useylm == 1) then 812 useylmgr1 = 1; optder = 1 813 end if 814 815 ABI_MALLOC(ylm_k, (npw_k, psps%mpsang**2 * psps%useylm)) 816 ABI_MALLOC(ylmgr1_k, (npw_k, 3+6*(optder/2), psps%mpsang**2*psps%useylm*useylmgr1)) 817 818 if (psps%useylm == 1) then 819 ! Fake MPI_type for sequential part. dummy_nband and nsppol1 are not used in sequential mode. 820 call initmpi_seq(mpienreg_seq) 821 dummy_nband = 0; npwarr = npw_k 822 call initylmg(cryst%gprimd, kg_k, kpoint, nkpt1, mpienreg_seq, psps%mpsang, npw_k, dummy_nband, nkpt1, & 823 npwarr, nsppol1, optder, cryst%rprimd, ylm_k, ylmgr1_k) 824 call destroy_mpi_enreg(mpienreg_seq) 825 end if 826 827 do idir=1,3 828 call self%htg(idir)%free() 829 830 ! Continue to initialize the Hamiltonian 831 call self%gs_hamkq(idir)%load_spin(spin, with_nonlocal=.true.) 832 call self%rf_hamkq(idir)%load_spin(spin, with_nonlocal=.true.) 833 834 ! We need ffnl1 and dkinpw for 3 dirs. Note that the Hamiltonian objects use pointers to keep a reference 835 ! to the output results of this routine. 836 ! This is the reason why we need to store the targets in self%htg 837 call getgh1c_setup(self%gs_hamkq(idir), self%rf_hamkq(idir), dtset, psps, kpoint, kpoint, idir, self%ipert, & ! In 838 cryst%natom, cryst%rmet, cryst%gprimd, cryst%gmet, istwf_k, npw_k, npw_k, & ! In 839 useylmgr1, kg_k, ylm_k, kg_k, ylm_k, ylmgr1_k, & ! In 840 self%htg(idir)%dkinpw, nkpg, nkpg1, self%htg(idir)%kpg_k, self%htg(idir)%kpg1_k, & ! Out 841 self%htg(idir)%kinpw1, self%htg(idir)%ffnlk, self%htg(idir)%ffnl1, & ! Out 842 self%htg(idir)%ph3d, self%htg(idir)%ph3d1) ! Out 843 end do 844 845 ABI_FREE(ylm_k) 846 ABI_FREE(ylmgr1_k) 847 848 end subroutine ddkop_setup_spin_kpoint
m_ddk/ddkop_t [ Types ]
NAME
ddkop_t
FUNCTION
This object provides a simplified interface to compute matrix elements of the velocity operator with the DFPT routines.
SOURCE
84 type,public :: ddkop_t 85 86 integer :: ipert 87 ! Perturbation type: natom + 1 88 89 integer :: inclvkb 90 ! Option for calculating the matrix elements of [Vnl,r]. 91 ! 0 to exclude commutator, 2 to include it 92 93 integer :: usepaw 94 ! 0 for NC, 1 for PAW 95 96 integer :: mpw 97 ! Maximum number of plane-waves over k-points (used to dimension arrays) 98 99 real(dp) :: kpoint(3) 100 ! K-point (set in setup_spin_kpoint) 101 102 real(dp) :: eig0nk 103 104 real(dp) :: dfpt_sciss = zero 105 106 real(dp) :: rprimd(3,3) 107 108 type(MPI_type),pointer :: mpi_enreg => null() 109 110 type(gs_hamiltonian_type) :: gs_hamkq(3) 111 112 type(rf_hamiltonian_type) :: rf_hamkq(3) 113 114 type(ham_targets_t), private :: htg(3) 115 ! Store arrays targetted by the hamiltonians. 116 117 real(dp), allocatable :: gh1c(:,:,:) 118 !gh1c, (2, mpw*nspinor, 3)) 119 120 real(dp), allocatable :: gs1c(:,:,:) 121 ! gs1c, (2, mpw*nspinor, 3*psps%usepaw)) 122 123 contains 124 125 procedure :: setup_spin_kpoint => ddkop_setup_spin_kpoint 126 ! Prepare application of dH/dk for given spin, k-point. 127 128 procedure :: apply => ddkop_apply 129 ! Apply dH/dk to input wavefunction. 130 131 procedure :: get_braket => ddkop_get_braket 132 ! Compute matrix element (complex results) in cartesian coords. 133 134 procedure :: get_vdiag => ddkop_get_vdiag 135 ! Compute diagonal matrix element (real) in cartesian coords. 136 137 procedure :: get_vnondiag => ddkop_get_vnondiag 138 ! Compute off diagonal matrix elements in cartesian coords. 139 140 procedure :: free => ddkop_free 141 ! Free memory. 142 143 end type ddkop_t 144 145 public :: ddkop_new ! Build object 146 147 !interface ddkop_t 148 ! procedure ddkop_new 149 !end interface ddkop_t
m_ddk/ddkstore_compute_ddk [ Functions ]
[ Top ] [ m_ddk ] [ Functions ]
NAME
ddkstore_compute_ddk
FUNCTION
Calculate the DDK matrix elements using the commutator formulation.
INPUTS
prefix: Prefix for output EVK file. Empty if output files are not wanted
SOURCE
211 subroutine ddkstore_compute_ddk(ds, wfk_path, prefix, dtset, psps, pawtab, ngfftc, comm) 212 213 !Arguments ------------------------------------ 214 !scalars 215 class(ddkstore_t),intent(inout) :: ds 216 character(len=*),intent(in) :: wfk_path, prefix 217 integer,intent(in) :: comm 218 type(dataset_type),intent(in) :: dtset 219 type(pseudopotential_type),intent(in) :: psps 220 type(pawtab_type),intent(in) :: pawtab(psps%ntypat*psps%usepaw) 221 !arrays 222 integer,intent(in) :: ngfftc(18) 223 224 !Local variables ------------------------------ 225 !scalars 226 integer,parameter :: master = 0 227 integer :: mband, nbcalc, nsppol, ib_v, ib_c, mpw, spin, nspinor, nkpt, nband_k, npw_k 228 integer :: ii, ik, bmin, bmax, istwf_k, idir, my_rank, nproc, ierr, bstop 229 real(dp) :: cpu, wall, gflops, cpu_all, wall_all, gflops_all 230 integer :: ncerr, ncid 231 character(len=500) :: msg 232 character(len=fnlen) :: fname 233 logical :: write_ncfile 234 type(wfd_t) :: wfd 235 type(vkbr_t) :: vkbr 236 type(ebands_t) :: ebands 237 type(crystal_t) :: cryst 238 type(hdr_type) :: tmp_hdr, hdr 239 type(ddkop_t) :: ddkop 240 type(wave_t),pointer :: wave_v, wave_c 241 !arrays 242 integer,allocatable :: distrib_mat(:,:,:,:), distrib_diago(:,:,:), nband(:,:), kg_k(:,:) 243 logical,allocatable :: bks_mask(:,:,:), keep_ur(:,:,:) 244 real(dp) :: kpt(3), vv(2, 3) 245 real(dp),allocatable :: cg_c(:,:), cg_v(:,:) 246 complex(dpc) :: vg(3), vr(3) 247 complex(gwpc),allocatable :: ihrc(:,:), ug_c(:), ug_v(:) 248 type(pawcprj_type),allocatable :: cwaveprj(:,:) 249 250 !************************************************************************ 251 252 my_rank = xmpi_comm_rank(comm); nproc = xmpi_comm_size(comm) 253 254 if (my_rank == master) call wrtout([std_out, ab_out], " Computation of velocity matrix elements (DDK)", newlines=1) 255 256 ABI_CHECK(psps%usepaw == 0, "PAW not implemented") 257 258 ! Get ebands and hdr from WFK file. 259 ebands = wfk_read_ebands(wfk_path, comm, out_hdr=hdr) 260 cryst = hdr%get_crystal() 261 262 ! Extract important dimensions from hdr% 263 nkpt = hdr%nkpt 264 nsppol = hdr%nsppol 265 nspinor = hdr%nspinor 266 mband = hdr%mband 267 268 ! Define band range 269 ! TODO: Perhaps one should allocate output arrays using bmin:bmax 270 ! and allow for nc output only if bmin == 1 and bmax == mband 271 if (ds%bmax == -1) ds%bmax = mband 272 273 if (ds%bmin < 1 .or. ds%bmin > mband .or. ds%bmax > mband .or. ds%bmin > ds%bmax) then 274 ABI_ERROR(sjoin("Invalid value for bmin, bmax", itoa(ds%bmin), itoa(ds%bmax), "with mband:", itoa(mband))) 275 end if 276 277 bmin = ds%bmin; bmax = ds%bmax 278 nbcalc = bmax - bmin + 1 279 write_ncfile = len_trim(prefix) > 0 280 if (write_ncfile .and. .not. (bmin == 1 .and. bmax == mband) ) then 281 write_ncfile = .False. 282 ABI_WARNING("Cannot write ncfile if .not. (bmin == 1 .and. bmax == mband)") 283 end if 284 285 if (my_rank == master) then 286 write(ab_out, "(a)")" Parameters extracted from the Abinit header:" 287 write(ab_out, "(a, f5.1)") ' ecut: ', hdr%ecut 288 write(ab_out, "(a, i0)") ' nkpt: ', nkpt 289 write(ab_out, "(a, i0)") ' mband: ', mband 290 write(ab_out, "(a, i0)") ' nsppol: ', nsppol 291 write(ab_out, "(a, i0)") ' nspinor: ', nspinor 292 write(ab_out, "(a, i0)") ' useylm: ', dtset%useylm 293 write(ab_out, "(a, i0)") ' inclvkb: ', dtset%inclvkb 294 write(ab_out, "(2(a, i0))")' bmin: ', bmin, ", bmax: ", bmax 295 if (ds%only_diago) then 296 write(ab_out, "(a)")' Computing diagonal matrix elements only' 297 else 298 write(ab_out, "(a)")' Computing diagonal and off-diagonal matrix elements' 299 end if 300 write(ab_out, "(2(a, i0))")' Between band index bmin: ', bmin, ", bmax: ", bmax 301 write(ab_out, "(a)")"" 302 end if 303 304 ! Create distribution of the wavefunctions mask. 305 ABI_MALLOC(nband, (nkpt, nsppol)) 306 ABI_MALLOC(keep_ur, (mband, nkpt, nsppol)) 307 ABI_MALLOC(bks_mask, (mband, nkpt, nsppol)) 308 keep_ur = .false.; bks_mask = .false.; nband = mband 309 310 if (ds%only_diago) then 311 ! Distribute k-points, spin and (b, b) diagonal over MPI processors. 312 ABI_MALLOC(distrib_diago, (bmin:bmax, nkpt, nsppol)) 313 distrib_diago = -1 314 315 ! Create bks_mask to load the wavefunctions. 316 ii = 0 317 do spin=1,nsppol 318 do ik=1,nkpt 319 do ib_v=bmin,bmax 320 ii = ii + 1; if (mod(ii, nproc) /= my_rank) cycle ! MPI parallelism. 321 distrib_diago(ib_v, ik, spin) = my_rank 322 bks_mask(ib_v, ik, spin) = .true. 323 end do 324 end do 325 end do 326 call wrtout(std_out, sjoin(" Rank: ", itoa(my_rank), "will treat", itoa(count(distrib_diago == my_rank)))) 327 328 else 329 ! Distribute k-points, spin and (b, b') pairs over the processors 330 ABI_MALLOC(distrib_mat, (bmin:bmax, bmin:bmax, nkpt, nsppol)) 331 call xmpi_distab(nproc, distrib_mat) 332 333 ! Create bks_mask to load the wavefunctions 334 do spin=1,nsppol 335 do ik=1,nkpt 336 ! Loop over v bands 337 do ib_v=bmin,bmax 338 ! Loop over c bands 339 do ib_c=bmin,bmax 340 if (distrib_mat(ib_c, ib_v, ik, spin) == my_rank) then 341 bks_mask(ib_v, ik, spin) = .true. 342 bks_mask(ib_c, ik, spin) = .true. 343 end if 344 end do 345 end do 346 end do 347 end do 348 349 call wrtout(std_out, sjoin(" Rank: ", itoa(my_rank), "will treat", itoa(count(distrib_mat == my_rank)))) 350 end if 351 352 ! Initialize distributed wavefunctions object 353 call wfd_init(wfd, cryst, pawtab, psps, keep_ur, mband, nband, nkpt, nsppol,& 354 bks_mask, dtset%nspden, nspinor, hdr%ecut, dtset%ecutsm, dtset%dilatmx, ebands%istwfk, ebands%kptns,& 355 ngfftc, dtset%nloalg, dtset%prtvol, dtset%pawprtvol, comm) 356 357 ABI_FREE(bks_mask) 358 ABI_FREE(keep_ur) 359 ABI_FREE(nband) 360 361 call wfd%print(header="Wavefunctions on the k-points grid") 362 363 ! Read wavefunctions from WFK file. 364 call wfd%read_wfk(wfk_path, iomode_from_fname(wfk_path)) 365 366 ! Allocate workspace arrays 367 mpw = maxval(wfd%npwarr) 368 ABI_MALLOC(kg_k, (3, mpw)) 369 ABI_MALLOC(ug_c, (mpw*nspinor)) 370 ABI_MALLOC(ug_v, (mpw*nspinor)) 371 if (dtset%useria /= 666) then 372 ABI_MALLOC(cg_c, (2, mpw*nspinor)) 373 ABI_MALLOC(cg_v, (2, mpw*nspinor)) 374 end if 375 376 ABI_MALLOC(cwaveprj, (0, 0)) 377 ABI_CALLOC(ds%dipoles, (3, 2, bmin:bmax, bmin:bmax, nkpt, nsppol)) 378 ABI_MALLOC(ihrc, (3, nspinor**2)) 379 380 if (ds%only_diago) then 381 ABI_CALLOC(ds%vdiago, (3, bmin:bmax, nkpt, nsppol)) 382 else 383 ABI_CALLOC(ds%vmat, (2, 3, bmin:bmax, bmin:bmax, nkpt, nsppol)) 384 end if 385 386 if (dtset%useria /= 666) then 387 ddkop = ddkop_new(dtset, cryst, pawtab, psps, wfd%mpi_enreg, mpw, wfd%ngfft) 388 !if (my_rank == master) call ddkop%print(ab_out) 389 end if 390 391 call cwtime(cpu_all, wall_all, gflops_all, "start") 392 393 do spin=1,nsppol 394 do ik=1,nkpt 395 396 ! Only do a subset a k-points 397 if (ds%only_diago) then 398 if (all(distrib_diago(:, ik, spin) /= my_rank)) cycle 399 else 400 if (all(distrib_mat(bmin:bmax, bmin:bmax, ik, spin) /= my_rank)) cycle 401 end if 402 call cwtime(cpu, wall, gflops, "start") 403 404 nband_k = wfd%nband(ik, spin) 405 istwf_k = wfd%istwfk(ik) 406 npw_k = wfd%npwarr(ik) 407 kpt = wfd%kibz(:,ik) 408 kg_k(:,1:npw_k) = wfd%kdata(ik)%kg_k 409 410 if (dtset%useria /= 666) then 411 call ddkop%setup_spin_kpoint(dtset, cryst, psps, spin, kpt, istwf_k, npw_k, kg_k) 412 else 413 ! Allocate KB form factors 414 ! Prepare term i <n,k|[Vnl,r]|n"k> 415 if (dtset%inclvkb /= 0) call vkbr_init(vkbr, cryst, psps, dtset%inclvkb, istwf_k, npw_k, kpt, kg_k) 416 end if 417 418 ! Loop over bands 419 do ib_v=bmin,bmax 420 if (ds%only_diago) then 421 if (distrib_diago(ib_v,ik,spin) /= my_rank) cycle 422 else 423 if (all(distrib_mat(:,ib_v,ik,spin) /= my_rank)) cycle 424 end if 425 426 if (dtset%useria /= 666) then 427 call wfd%copy_cg(ib_v, ik, spin, cg_v) 428 call ddkop%apply(ebands%eig(ib_v, ik, spin), npw_k, wfd%nspinor, cg_v, cwaveprj) 429 else 430 ABI_CHECK(wfd%get_wave_ptr(ib_v, ik, spin, wave_v, msg) == 0, msg) 431 ug_v(1:npw_k*nspinor) = wave_v%ug 432 end if 433 434 ! Loop over bands 435 bstop = bmax; if (ds%only_diago) bstop = ib_v 436 do ib_c=ib_v,bstop 437 if (.not. ds%only_diago) then 438 if (distrib_mat(ib_c, ib_v, ik, spin) /= my_rank) cycle 439 end if 440 441 if (dtset%useria /= 666) then 442 call wfd%copy_cg(ib_c, ik, spin, cg_c) 443 vv = ddkop%get_braket(ebands%eig(ib_c, ik, spin), istwf_k, npw_k, nspinor, cg_c, mode=ds%mode) 444 !if (ib_v == ib_c) vv(2, :) = zero 445 446 if (ds%only_diago) then 447 ds%vdiago(:,ib_c,ik,spin) = vv(1, :) 448 else 449 ds%vmat(:,:,ib_c,ib_v,ik,spin) = vv 450 ! Hermitian conjugate 451 if (ib_v /= ib_c) then 452 ds%vmat(1,:,ib_v,ib_c,ik,spin) = vv(1, :) 453 ds%vmat(2,:,ib_v,ib_c,ik,spin) = -vv(2, :) 454 end if 455 end if 456 457 do idir=1,3 458 ds%dipoles(idir,:,ib_c,ib_v,ik,spin) = vv(:, idir) 459 ! Hermitian conjugate 460 if (ib_v /= ib_c) ds%dipoles(idir,:,ib_v,ib_c,ik,spin) = [vv(1, idir), -vv(2, idir)] 461 end do 462 463 else 464 ABI_CHECK(wfd%get_wave_ptr(ib_c, ik, spin, wave_c, msg) == 0, msg) 465 ug_c(1:npw_k*nspinor) = wave_c%ug 466 467 ! Calculate matrix elements of i[H,r] for NC pseudopotentials. 468 ihrc = nc_ihr_comm(vkbr, cryst, psps, npw_k, nspinor, istwf_k, dtset%inclvkb, kpt, ug_c, ug_v, kg_k) 469 470 ! HM: 24/07/2018 471 ! Transform dipoles to be consistent with results from DFPT 472 ! Perturbations with DFPT are along the reciprocal lattice vectors 473 ! Perturbations with Commutator are along real space lattice vectors 474 ! dot(A, DFPT) = X 475 ! dot(B, COMM) = X 476 ! B = 2 pi (A^{-1})^T => 477 ! dot(B^T B,COMM) = 2 pi DFPT 478 vr = (2*pi)*(2*pi)*sum(ihrc(:,:),dim=2) 479 vg(1) = dot_product(Cryst%gmet(1,:), vr) 480 vg(2) = dot_product(Cryst%gmet(2,:), vr) 481 vg(3) = dot_product(Cryst%gmet(3,:), vr) 482 483 ! Save matrix elements of i*r in the IBZ 484 ds%dipoles(:,1,ib_c,ib_v,ik,spin) = real(vg, kind=dp) 485 ds%dipoles(:,1,ib_v,ib_c,ik,spin) = real(vg, kind=dp) ! Hermitian conjugate 486 if (ib_v == ib_c) then 487 ds%dipoles(:,2,ib_c,ib_v,ik,spin) = zero 488 ds%dipoles(:,2,ib_v,ib_c,ik,spin) = zero 489 else 490 ds%dipoles(:,2,ib_c,ib_v,ik,spin) = aimag(vg) 491 ds%dipoles(:,2,ib_v,ib_c,ik,spin) = -aimag(vg) ! Hermitian conjugate 492 end if 493 end if 494 495 end do 496 end do 497 498 ! Free KB form factors 499 call vkbr_free(vkbr) 500 501 if (nkpt < 1000 .or. (nkpt > 1000 .and. mod(ik, 200) == 0) .or. ik <= nproc) then 502 write(msg,'(2(a,i0),a)')" k-point [", ik, "/", nkpt, "]" 503 call cwtime_report(msg, cpu, wall, gflops) 504 end if 505 506 end do ! k-points 507 end do ! spin 508 509 call cwtime_report(msg, cpu_all, wall_all, gflops_all) 510 511 ABI_FREE(ug_c) 512 ABI_FREE(ug_v) 513 ABI_FREE(kg_k) 514 ABI_FREE(ihrc) 515 ABI_FREE(cwaveprj) 516 ABI_SFREE(distrib_mat) 517 ABI_SFREE(distrib_diago) 518 519 if (dtset%useria /= 666) then 520 ABI_FREE(cg_c) 521 ABI_FREE(cg_v) 522 call ddkop%free() 523 end if 524 525 ! Gather the k-points computed by all processes 526 call xmpi_sum_master(ds%dipoles, master, comm, ierr) 527 528 if (ds%only_diago) then 529 call xmpi_sum_master(ds%vdiago, master, comm, ierr) 530 else 531 call xmpi_sum_master(ds%vmat, master, comm, ierr) 532 end if 533 534 ! Write matrix elements to disk. 535 536 ! Output EVK file in netcdf format. 537 if (my_rank == master .and. write_ncfile) then 538 ! Have to build hdr on k-grid with info about perturbation. 539 call hdr_copy(hdr, tmp_hdr) 540 tmp_hdr%qptn = zero 541 542 !fname = strcat(prefix, "NEW_EVK.nc") 543 !call wrtout(ab_out, sjoin("- Writing file: ", fname)) 544 !NCF_CHECK_MSG(nctk_open_create(ncid, fname, xmpi_comm_self), "Creating EVK.nc file") 545 !tmp_hdr%pertcase = 0 546 !NCF_CHECK(tmp_hdr%ncwrite(ncid, 43, nc_define=.True.)) 547 !NCF_CHECK(cryst%ncwrite(ncid)) 548 !NCF_CHECK(ebands_ncwrite(ebands, ncid)) 549 !if (ds%only_diago) then 550 ! ncerr = nctk_def_arrays(ncid, [ & 551 ! nctkarr_t('vred_diagonal', "dp", "three, max_number_of_states, number_of_kpoints, number_of_spins")], defmode=.True.) 552 !else 553 ! ncerr = nctk_def_arrays(ncid, [ nctkarr_t('vred_matrix', "dp", & 554 ! "two, three, max_number_of_states, max_number_of_states, number_of_kpoints, number_of_spins")], defmode=.True.) 555 !end if 556 !NCF_CHECK(ncerr) 557 !NCF_CHECK(nctk_set_datamode(ncid)) 558 !if (ds%only_diago) then 559 ! NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "vred_diagonal"), ds%vdiago)) 560 !else 561 ! NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "vred_matrix"), ds%vmat)) 562 !end if 563 !NCF_CHECK(nf90_close(ncid)) 564 565 do ii=1,3 566 fname = strcat(prefix, '_', itoa(ii), "_EVK.nc") 567 call wrtout(ab_out, sjoin("- Writing EVK file: ", fname, "for reduced direction:", itoa(ii))) 568 NCF_CHECK_MSG(nctk_open_create(ncid, fname, xmpi_comm_self), "Creating EVK.nc file") 569 tmp_hdr%pertcase = 3 * cryst%natom + ii 570 NCF_CHECK(tmp_hdr%ncwrite(ncid, 43, nc_define=.True.)) 571 NCF_CHECK(cryst%ncwrite(ncid)) 572 NCF_CHECK(ebands_ncwrite(ebands, ncid)) 573 ncerr = nctk_def_arrays(ncid, [ & 574 nctkarr_t('h1_matrix_elements', "dp", & 575 "two, max_number_of_states, max_number_of_states, number_of_kpoints, number_of_spins")], defmode=.True.) 576 NCF_CHECK(ncerr) 577 NCF_CHECK(nctk_set_datamode(ncid)) 578 NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "h1_matrix_elements"), ds%dipoles(ii,:,:,:,:,:))) 579 NCF_CHECK(nf90_close(ncid)) 580 end do 581 call tmp_hdr%free() 582 end if 583 584 if (my_rank == master .and. dtset%prtvol > 0) then 585 write(ab_out, "(2a)")ch10,"Writing velocity matrix elements (only diagonal terms, real part) for testing purpose:" 586 do spin=1,nsppol 587 do ik=1,min(nkpt, 4) 588 write(ab_out, "(2(a,1x,i0),2x,2a)")"For spin: ", spin, ", ikbz: ", ik, ", kpt: ", trim(ktoa(wfd%kibz(:,ik))) 589 do ib_c=bmin,min(bmin+8, bmax) 590 write(ab_out, "(3(es16.6,2x))") ds%dipoles(:,1,ib_c,ib_c,ik,spin) 591 end do 592 write(ab_out,*)"" 593 !do ib_c=bmin,min(bmin+8, bmax) 594 ! write(ab_out, "(a, 6(es16.6,2x))")"Sum_k: ", sum(ds%dipoles(:,:,ib_c,ib_c,:,spin), dim=3) / nkpt 595 !end do 596 end do 597 end do 598 end if 599 600 ! Free memory 601 call wfd%free() 602 call ebands_free(ebands) 603 call cryst%free() 604 call hdr%free() 605 606 ! Block all procs here so that we know output files are available when code returns. 607 call xmpi_barrier(comm) 608 609 end subroutine ddkstore_compute_ddk
m_ddk/ddkstore_free [ Functions ]
[ Top ] [ m_ddk ] [ Functions ]
NAME
ddkstore_free
FUNCTION
Free memory
INPUTS
SOURCE
625 subroutine ddkstore_free(self) 626 627 !Arguments ------------------------------------ 628 !scalars 629 class(ddkstore_t),intent(inout) :: self 630 631 !************************************************************************ 632 633 ABI_SFREE(self%vdiago) 634 ABI_SFREE(self%vmat) 635 ABI_SFREE(self%dipoles) 636 637 end subroutine ddkstore_free
m_ddk/ddkstore_t [ Types ]
NAME
ddkstore_t
FUNCTION
This object stores the matrix elements of the velocity operator computed with the DFPT routines.
SOURCE
161 type,public :: ddkstore_t 162 163 integer :: bmin = 1, bmax = -1 164 ! Min and max band index 165 166 character(len=50) :: mode = "reduced" 167 ! "cart" or "reduced" 168 169 logical :: only_diago = .False. 170 ! True if we are computing only the diagonal elements 171 172 real(dp),allocatable :: dipoles(:,:,:,:,:,:) 173 ! (3, 2, mband, mband, nkpt, nsppol)) 174 175 real(dp),allocatable :: vdiago(:,:,:,:) 176 ! (3, bmin:bmax, nkpt, nsppol) 177 178 real(dp),allocatable :: vmat(:,:,:,:,:,:) 179 ! (2, 3, bmin:bmax, bmin:bmax, nkpt, nsppol)) 180 181 contains 182 183 procedure :: compute_ddk => ddkstore_compute_ddk 184 ! Calculate DDK matrix elements (diago or full b,b' matrix). 185 ! Return results in datatype. Optionally, save results to disk in EVK format. 186 187 procedure :: free => ddkstore_free 188 ! Free memory. 189 190 end type ddkstore_t
m_ddk/ham_targets_free [ Functions ]
[ Top ] [ m_ddk ] [ Functions ]
NAME
FUNCTION
INPUTS
SOURCE
1140 subroutine ham_targets_free(self) 1141 1142 !Arguments ------------------------------------ 1143 !scalars 1144 class(ham_targets_t),intent(inout) :: self 1145 1146 !************************************************************************ 1147 1148 ABI_SFREE(self%ffnlk) 1149 ABI_SFREE(self%ffnl1) 1150 ABI_SFREE(self%kpg_k) 1151 ABI_SFREE(self%kpg1_k) 1152 ABI_SFREE(self%dkinpw) 1153 ABI_SFREE(self%kinpw1) 1154 ABI_SFREE(self%ph3d) 1155 ABI_SFREE(self%ph3d1) 1156 1157 end subroutine ham_targets_free