TABLE OF CONTENTS
ABINIT/interpolate_gkk [ Functions ]
NAME
interpolate_gkk
FUNCTION
This routine interpolates the gkk matrices for all q vectors between points on the full kpt_phon grid.
COPYRIGHT
Copyright (C) 2004-2018 ABINIT group (MVer) This file is distributed under the terms of the GNU General Public Licence, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt . For the initials of contributors, see ~abinit/doc/developers/contributors.txt .
INPUTS
elph_ds = elphon datastructure with data and dimensions kpt_phon = coordinates of all kpoints close to the FS
OUTPUT
elph_ds = modified gkq
NOTES
inspired to some extent by epcouple.f from the DecAFT package by J. Kay Dewhurst most inputs taken from mkifc.f in anaddb set ifcflag 1 such that the IFC are calculated in atmfrc prior to calling elphon
PARENTS
get_all_gkk2
CHILDREN
ftgkk,ifc_fourq,wrap2_pmhalf,zhpev
SOURCE
38 #if defined HAVE_CONFIG_H 39 #include "config.h" 40 #endif 41 42 #include "abi_common.h" 43 44 45 subroutine interpolate_gkk(crystal,ifc,elph_ds,kpt_phon) 46 47 use defs_basis 48 use defs_elphon 49 use m_profiling_abi 50 use m_errors 51 52 use m_numeric_tools, only : wrap2_pmhalf 53 use m_crystal, only : crystal_t 54 use m_ifc, only : ifc_type, ifc_fourq 55 56 !This section has been created automatically by the script Abilint (TD). 57 !Do not modify the following lines by hand. 58 #undef ABI_FUNC 59 #define ABI_FUNC 'interpolate_gkk' 60 use interfaces_77_ddb, except_this_one => interpolate_gkk 61 !End of the abilint section 62 63 implicit none 64 65 !Arguments ------------------------------------ 66 !scalars 67 type(crystal_t),intent(in) :: crystal 68 type(ifc_type),intent(in) :: ifc 69 type(elph_type),intent(inout) :: elph_ds 70 !arrays 71 real(dp),intent(in) :: kpt_phon(3,elph_ds%k_phon%nkpt) 72 73 !Local variables------------------------------- 74 ! output variables for dfpt_phfrq 75 ! variables for zhpev 76 ! variables for phonon interpolation 77 !scalars 78 integer :: i1,i2,ikpt_phon2,iFSqpt,ib1,ib2,ier,ii 79 integer :: iost,isppol,qtor,natom 80 integer :: sz1,sz2,sz3,sz4,unit_gkkp 81 real(dp) :: qphnrm,res 82 !character(len=500) :: msg 83 !arrays 84 real(dp) :: gprim(3,3) 85 real(dp) :: displ(2,elph_ds%nbranch,elph_ds%nbranch),eigval(3*crystal%natom) 86 real(dp) :: eigvec(3*3*crystal%natom*3*crystal%natom) 87 real(dp) :: pheigvec(2*elph_ds%nbranch*elph_ds%nbranch) 88 real(dp) :: phfrq_tmp(elph_ds%nbranch),qphon(3),redkpt(3) 89 real(dp),allocatable :: gkk2_diag_tmp(:,:,:,:),gkk2_tmp(:,:,:,:,:,:,:) 90 real(dp),allocatable :: matrx(:,:),zhpev1(:,:) 91 real(dp),allocatable :: zhpev2(:) 92 93 ! ************************************************************************* 94 95 ! 96 !NOTE: mjv 18/5/2008 reverted to old style of ftgkk with all kpt done together. 97 !may want to modify this later to use the new cleaner format with 1 FT at a 98 !time. 99 ! 100 write(std_out,*) 'interpolate_gkk : enter' 101 102 natom = crystal%natom 103 gprim = ifc%gprim 104 105 if (elph_ds%nsppol /= 1) then 106 MSG_ERROR("interpolate_gkk not coded with nsppol>1 yet") 107 end if 108 isppol = 1 109 110 111 !------------------------------------------------------ 112 !complete dynamical matrices for all qpts between points 113 !on full kpt grid (interpolation from IFC) 114 !------------------------------------------------------ 115 116 sz1=elph_ds%ngkkband;sz2=elph_ds%nbranch 117 sz3=elph_ds%k_phon%nkpt;sz4=elph_ds%nFSband 118 !allocate (gkk_tmp(2,sz1,sz1,sz2,sz2,1,1)) 119 !DEBUG 120 !allocate (gkk_tmp_full(2,sz1,sz1,sz2,elph_ds%nFSband,sz3)) 121 !allocate (gkk_tmp_full(2,s2,sz4,sz4,sz3)) 122 !ENDDEBUG 123 ABI_ALLOCATE(gkk2_tmp,(2,sz1,sz1,sz2,sz2,sz3,1)) 124 ABI_ALLOCATE(gkk2_diag_tmp,(sz1,sz1,sz2,sz3)) 125 ABI_ALLOCATE(zhpev1,(2,2*3*natom-1)) 126 ABI_ALLOCATE(zhpev2,(3*3*natom-2)) 127 ABI_ALLOCATE(matrx,(2,(3*natom*(3*natom+1))/2)) 128 129 qphnrm = one 130 !in this part use the inverse Fourier transform to get 1 (arbitrary) qpt at a 131 !time 132 ii = 0 133 qtor = 0 134 unit_gkkp = 150 135 open (unit=unit_gkkp,file='gkkp_file_ascii',form='formatted',status='unknown',iostat=iost) 136 if (iost /= 0) then 137 MSG_ERROR("error opening gkkpfile as new") 138 end if 139 140 !loop over all FS pairs. 141 !do ikpt1=1,elph_ds%k_phon%nkptirr 142 !do iFSqpt=1,elph_ds%k_phon%nkpt 143 144 ! 145 !this should run through the sparse mesh of 2x2x2 kpoints 146 ! 147 do iFSqpt=1,elph_ds%k_phon%nkpt 148 res = 2.0_dp*(kpt_phon(1,iFSqpt)+one) 149 if (abs(res-int(res)) > tol10) cycle 150 res = 2.0_dp*(kpt_phon(2,iFSqpt)+one) 151 if (abs(res-int(res)) > tol10) cycle 152 res = 2.0_dp*(kpt_phon(3,iFSqpt)+one) 153 if (abs(res-int(res)) > tol10) cycle 154 155 ! do ikpt1=1,1 156 ! 157 ! NOTE: should be very easy to parallelize! 158 ! 159 ! write(std_out,*) ' interpolate_gkk : ikpt1 = ',ikpt1, ' / ', elph_ds%k_phon%nkptirr 160 write(std_out,*) ' interpolate_gkk : ikpt1 = ',iFSqpt, ' / ', elph_ds%k_phon%nkpt 161 162 ! DEBUG 163 ! write(std_out,*) ' interpolate_gkk : Warning debug version' 164 ! cycle 165 ! ENDDEBUG 166 167 gkk2_tmp(:,:,:,:,:,:,:) = zero 168 169 ! qphon = 1 - 2 ie. 1 = 2+qphon 170 qphon(:) = kpt_phon(:,iFSqpt) 171 172 ! shouldnt be necessary here, but oh well 173 call wrap2_pmhalf(qphon(1),redkpt(1),res) 174 call wrap2_pmhalf(qphon(2),redkpt(2),res) 175 call wrap2_pmhalf(qphon(3),redkpt(3),res) 176 177 qphon(:) = redkpt(:) 178 redkpt(1) = qphon(1)*gprim(1,1)+qphon(2)*gprim(1,2)+qphon(3)*gprim(1,3) 179 redkpt(2) = qphon(1)*gprim(2,1)+qphon(2)*gprim(2,2)+qphon(3)*gprim(2,3) 180 redkpt(3) = qphon(1)*gprim(3,1)+qphon(2)*gprim(3,2)+qphon(3)*gprim(3,3) 181 write (unit_gkkp,*) 'qp= ', redkpt 182 183 call ifc_fourq(ifc,crystal,qphon,phfrq_tmp,displ,out_eigvec=pheigvec) 184 write (unit_gkkp,*) phfrq_tmp(:)*Ha_cmm1 185 186 ii = ii+1 187 ! if(ii > 0 .and. ii < 1000) write(std_out,'(a,i5,3E16.6,2x)') & 188 ! & ' wrote phfrq_tmp for time ', ii, phfrq_tmp 189 ! end if 190 191 ! phonon eigenvectors are in eigvec 192 ! real and imaginary parts 193 ! phonon displacements = eigvec/sqrt(M_i) are in displ 194 ! real and imaginary parts 195 196 ! DEBUG 197 ! test: uniform phonon frequency 198 ! phfrq_tmp(:) = 0.0001_dp 199 ! ENDDEBUG 200 201 ! FT gamma matrices for all kpt_phon points, and 202 ! for qpoint = qphon(:) = kpt_phon(ikpt_phon) 203 204 call ftgkk(ifc%wghatm,gkk2_tmp,elph_ds%gkk_rpt,elph_ds%gkqwrite,& 205 & elph_ds%gkk_rptwrite,gprim,1,& 206 & natom,elph_ds%k_phon%nkpt,elph_ds%ngkkband,elph_ds%k_phon%nkpt,1,ifc%nrpt,elph_ds%nsppol,& 207 & qtor,ifc%rpt,qphon,elph_ds%unit_gkk_rpt,elph_ds%unitgkq) 208 209 ! NOTE: Normally the eigenvectors of the gkk2_tmp should be the same as eigvec 210 211 ! Diagonalize gamma matrices at qpoint (complex matrix) for all kpt_phon. 212 ! Copied from dfpt_phfrq 213 do ikpt_phon2=1,elph_ds%k_phon%nkpt 214 res = 8.0_dp*(kpt_phon(1,ikpt_phon2)+one) 215 if (abs(res-int(res)) > tol10) cycle 216 res = 8.0_dp*(kpt_phon(2,ikpt_phon2)+one) 217 if (abs(res-int(res)) > tol10) cycle 218 res = 8.0_dp*(kpt_phon(3,ikpt_phon2)+one) 219 if (abs(res-int(res)) > tol10) cycle 220 221 write (unit_gkkp,*) 'kp= ', kpt_phon(:,ikpt_phon2) 222 223 do ib1=1,elph_ds%ngkkband 224 do ib2=1,elph_ds%ngkkband 225 ier=0 226 ii=1 227 do i2=1,3*natom 228 do i1=1,i2 229 matrx(1,ii)=gkk2_tmp(1,ib1,ib2,i1,i2,ikpt_phon2,1) 230 matrx(2,ii)=gkk2_tmp(2,ib1,ib2,i1,i2,ikpt_phon2,1) 231 ii=ii+1 232 end do 233 end do 234 call ZHPEV ('N','U',3*natom,matrx,eigval,eigvec,3*natom,zhpev1,& 235 & zhpev2,ier) 236 237 gkk2_diag_tmp(ib2,ib1,:,ikpt_phon2) = eigval(:) 238 do i1=1,3*natom 239 write (unit_gkkp,*) elph_ds%minFSband-1+ib1,elph_ds%minFSband-1+ib2,i1,& 240 & eigval(i1) 241 end do 242 end do 243 end do 244 end do 245 246 if (elph_ds%gkk2write == 1) then 247 write(std_out,*) 'WARNING COMMENTED WRITE TO BINARY FILE!!!' 248 ! write (elph_ds%unit_gkk2,REC=iFSqpt) gkk2_diag_tmp(:,:,:,:) 249 write(std_out,'(a,i4,4(2E16.6,2x))') ' gkk2 loop ', & 250 & iFSqpt,gkk2_diag_tmp(1,1,:,1:2),gkk2_diag_tmp(1,1,:,elph_ds%k_phon%nkpt-1:elph_ds%k_phon%nkpt) 251 ! & ikpt1,gkk2_tmp(:,1,1,1,1,1:2),gkk2_tmp(:,1,1,elph_ds%k_phon%nkpt-1:elph_ds%k_phon%nkpt) 252 else if (elph_ds%gkk2write == 0) then 253 elph_ds%gkk2(:,:,:,:,iFSqpt,isppol) = gkk2_diag_tmp(:,:,:,:) 254 ! elph_ds%gkk2(:,:,:,:,ikpt1) = gkk2_tmp 255 write(std_out,*) ' interpolate_gkk : gkk2(b=1,b=1,:,kpt=1,iFSqpt) = ' 256 write(std_out,*) gkk2_diag_tmp(1,1,:,1) 257 end if 258 259 end do 260 !end do on iFSqpt 261 262 ABI_DEALLOCATE(matrx) 263 ABI_DEALLOCATE(zhpev1) 264 ABI_DEALLOCATE(zhpev2) 265 266 end subroutine interpolate_gkk