TABLE OF CONTENTS
ABINIT/complete_gamma_tr [ Functions ]
NAME
complete_gamma_tr
FUNCTION
Use the set of special q points calculated by the Monkhorst & Pack Technique. Check if all the informations for the q points are present in the input gamma transport matrices. Generate the gamma transport matrices (already summed over the FS) of the set of q points which samples homogeneously the entire Brillouin zone.
COPYRIGHT
Copyright (C) 2009-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
crystal<crystal_t>=data type gathering info on the crystalline structure. ep_scalprod= flag for scalar product of gkk with phonon displacement vectors nbranch=number of phonon branches = 3*natom nqptirred=nqpt irred BZ nqpt_full=nqpt full BZ nsppol=number of spins qirredtofull= mapping irred to full qpoints qpttoqpt = qpoint index mapping under symops
OUTPUT
gamma_qpt_tr = in/out: set of gamma matrix elements completed and symmetrized gamma_qpt_tr(2,9,nbranch*nbranch,nsppol,nqpt_full)
PARENTS
elphon
CHILDREN
dgemm
SOURCE
43 #if defined HAVE_CONFIG_H 44 #include "config.h" 45 #endif 46 47 #include "abi_common.h" 48 49 subroutine complete_gamma_tr(crystal,ep_scalprod,nbranch,nqptirred,nqpt_full,nsppol,gamma_qpt_tr,qirredtofull,qpttoqpt) 50 51 use defs_basis 52 use defs_elphon 53 use m_profiling_abi 54 use m_linalg_interfaces 55 use m_errors 56 57 use m_crystal, only : crystal_t 58 59 !This section has been created automatically by the script Abilint (TD). 60 !Do not modify the following lines by hand. 61 #undef ABI_FUNC 62 #define ABI_FUNC 'complete_gamma_tr' 63 !End of the abilint section 64 65 implicit none 66 67 !Arguments ------------------------------------ 68 !scalars 69 integer, intent(in) :: nbranch,nqptirred,nqpt_full,nsppol, ep_scalprod 70 type(crystal_t),intent(in) :: crystal 71 !arrays 72 integer,intent(in) :: qpttoqpt(2,crystal%nsym,nqpt_full) 73 integer,intent(in) :: qirredtofull(nqptirred) 74 real(dp), intent(inout) :: gamma_qpt_tr(2,9,nbranch*nbranch,nsppol,nqpt_full) 75 76 !Local variables------------------------------- 77 !scalars 78 integer :: ieqqpt,ii,iqpt,isppol,isym 79 integer :: itim,jj,kk,ll,neqqpt 80 integer :: iatom,ancestor_iatom 81 integer :: iqpt_fullbz,imode, itensor 82 real(dp),parameter :: tol=2.d-8 83 !arrays 84 integer :: symrel(3,3,crystal%nsym),symrec(3,3,crystal%nsym) 85 integer :: symmetrized_qpt(nqpt_full) 86 integer :: gkk_flag(nbranch,nbranch,nsppol,nqpt_full) 87 real(dp) :: gprimd(3,3),rprimd(3,3) 88 real(dp) :: ss(3,3), sscart(3,3) 89 real(dp) :: tmp_mat(2,nbranch,nbranch) 90 real(dp) :: tmp_mat2(2,nbranch,nbranch) 91 real(dp) :: tmp_tensor(2,3,3) 92 real(dp) :: tmp_tensor2(2,3,3) 93 real(dp) :: ss_allatoms(nbranch,nbranch) 94 real(dp) :: c_one(2), c_zero(2) 95 real(dp),allocatable :: gkk_qpt_new(:,:,:,:),gkk_qpt_tmp(:,:,:,:) 96 97 ! ********************************************************************* 98 99 c_one = (/one,zero/) 100 c_zero = (/zero,zero/) 101 102 gprimd = crystal%gprimd 103 rprimd = crystal%rprimd 104 105 symrec = crystal%symrec 106 symrel = crystal%symrel 107 108 !Generation of the gkk matrices relative to the q points 109 !of the set which samples the entire Brillouin zone 110 111 !set up flags for gamma_qpt matrices we have 112 gkk_flag = -1 113 do iqpt=1,nqptirred 114 iqpt_fullbz = qirredtofull(iqpt) 115 gkk_flag(:,:,:,iqpt_fullbz) = 1 116 end do 117 118 symmetrized_qpt(:) = -1 119 ! isppol=1 120 121 ABI_ALLOCATE(gkk_qpt_new,(2,9,nbranch*nbranch, nsppol)) 122 ABI_ALLOCATE(gkk_qpt_tmp,(2,9,nbranch*nbranch, nsppol)) 123 124 do iqpt=1,nqpt_full 125 126 ! Already symmetrized? 127 if (symmetrized_qpt(iqpt) == 1) cycle 128 129 gkk_qpt_new(:,:,:,:) = zero 130 131 ! loop over qpoints equivalent to iqpt 132 neqqpt=0 133 ! do not use time reversal symmetry to complete the qpoints: 134 ! do not know what happens to the gamma matrices 135 136 do itim=1,2 137 do isym=1,crystal%nsym 138 ! ieqqpt is sent onto iqpt by itim/isym 139 ieqqpt = qpttoqpt(itim,isym,iqpt) 140 141 if (gkk_flag(1,1,1,ieqqpt) == -1) cycle 142 ! if we have information on this qpt 143 ! iqpt is equivalent to ieqqpt: get it from file or memory 144 gkk_qpt_tmp(:,:,:,:) = gamma_qpt_tr(:,:,:,:,ieqqpt) 145 146 neqqpt=neqqpt+1 147 148 ! 149 ! MJV note 02/2010: 150 ! the correspondence of symrel and symrec in the different cases, symmetrizing there 151 ! and back, has been fixed in the cases with and without scalprod (ie cartesian 152 ! and reduced real space coordinates) with respect to a calculation with no symmetries 153 ! I believe everything is settled, but still do not know why the 2 versions of the ss 154 ! matrices here use different rel/rec, instead of just being multiplied by the rprim gprim... 155 ! 156 do ii=1,3 157 do jj=1,3 158 sscart(ii,jj)=0.0_dp 159 do kk=1,3 160 do ll=1,3 161 sscart(ii,jj)=sscart(ii,jj)+rprimd(ii,kk)*symrel(kk,ll,isym)*gprimd(ll,jj) 162 ! sscart(ii,jj)=sscart(ii,jj)+rprimd(ii,kk)*symrel(kk,ll,isym)*gprimd(ll,jj) 163 end do 164 end do 165 end do 166 end do 167 if (ep_scalprod==1) then 168 do ii=1,3 169 do jj=1,3 170 ss(ii,jj)=0.0_dp 171 do kk=1,3 172 do ll=1,3 173 ss(ii,jj)=ss(ii,jj)+rprimd(ii,kk)*symrel(kk,ll,isym)*gprimd(ll,jj) 174 end do 175 end do 176 end do 177 end do 178 else 179 do ii=1,3 180 do jj=1,3 181 ss(ii,jj) = symrec(ii,jj,isym) 182 end do 183 end do 184 end if 185 186 ss_allatoms(:,:) = zero 187 do iatom=1,crystal%natom 188 ancestor_iatom = crystal%indsym(4,isym,iatom) 189 ss_allatoms((ancestor_iatom-1)*3+1:(ancestor_iatom-1)*3+3,& 190 & (iatom-1)*3+1: (iatom-1)*3+3) = ss(1:3,1:3) 191 end do 192 193 194 ! NOTE ssinv(ii,jj)=ssinv(ii,jj)+gprimd(ii,kk)*rprimd(jj,ll)*symrec(ll,kk,isym) 195 196 do isppol=1,nsppol 197 198 ! for each tensor component, rotate the cartesian directions of phonon modes 199 do itensor = 1, 9 200 ! multiply by the ss matrices 201 tmp_mat2(:,:,:) = zero 202 tmp_mat(:,:,:) = reshape(gkk_qpt_tmp(:,itensor,:,isppol),& 203 & (/2,nbranch,nbranch/)) 204 call DGEMM ('N','N',nbranch,nbranch,nbranch,& 205 & one,ss_allatoms,nbranch,tmp_mat(1,:,:),nbranch,zero,& 206 & tmp_mat2(1,:,:),nbranch) 207 call DGEMM ('N','N',nbranch,nbranch,nbranch,& 208 & one,ss_allatoms,nbranch,tmp_mat(2,:,:),nbranch,zero,& 209 & tmp_mat2(2,:,:),nbranch) 210 211 call DGEMM ('N','T',nbranch,nbranch,nbranch,& 212 & one,tmp_mat2(1,:,:),nbranch,ss_allatoms,nbranch,zero,& 213 & tmp_mat(1,:,:),nbranch) 214 call DGEMM ('N','T',nbranch,nbranch,nbranch,& 215 & one,tmp_mat2(2,:,:),nbranch,ss_allatoms,nbranch,zero,& 216 & tmp_mat(2,:,:),nbranch) 217 218 gkk_qpt_tmp(:,itensor,:,isppol) = reshape (tmp_mat, (/2,nbranch*nbranch/)) 219 end do ! itensor 220 221 ! for each cartesian direction/phonon mode, rotate the tensor components 222 do imode = 1, nbranch*nbranch 223 tmp_tensor2(:,:,:) = zero 224 tmp_tensor(:,:,:) = reshape(gkk_qpt_tmp(:,:,imode,isppol),& 225 & (/2,3,3/)) 226 call DGEMM ('N','N',3,3,3,& 227 & one,sscart,3,tmp_tensor(1,:,:),3,zero,& 228 & tmp_tensor2(1,:,:),3) 229 call DGEMM ('N','T',3,3,3,& 230 & one,tmp_tensor2(1,:,:),3,sscart,3,zero,& 231 & tmp_tensor(1,:,:),3) 232 233 call DGEMM ('N','N',3,3,3,& 234 & one,sscart,3,tmp_tensor(2,:,:),3,zero,& 235 & tmp_tensor2(2,:,:),3) 236 call DGEMM ('N','T',3,3,3,& 237 & one,tmp_tensor2(2,:,:),3,sscart,3,zero,& 238 & tmp_tensor(2,:,:),3) 239 240 gkk_qpt_tmp(:,:,imode,isppol) = reshape (tmp_tensor, (/2,9/)) ! modified by BX 241 end do ! imode 242 243 ! add to gkk_qpt_new 244 gkk_qpt_new(:,:,:,isppol) = gkk_qpt_new(:,:,:,isppol) + gkk_qpt_tmp(:,:,:,isppol) 245 246 end do ! end isppol do 247 248 end do ! end isym do 249 end do ! end itim do 250 251 ABI_CHECK(neqqpt>0,'no q-points found equivalent to iqpt ') 252 253 ! divide by number of equivalent qpts found 254 gkk_qpt_new = gkk_qpt_new/neqqpt 255 256 257 ! copy the symmetrized version into all the equivalent qpoints, appropriately transformed 258 do itim=1,2 259 do isym=1,crystal%nsym 260 ! ieqqpt is sent onto iqpt by itim/isym 261 ieqqpt = qpttoqpt(itim,isym,iqpt) 262 263 if (symmetrized_qpt(ieqqpt) /= -1) cycle 264 gkk_qpt_tmp = zero 265 266 ! use symrec matrices to get inverse transform from isym^{-1} 267 do ii=1,3 268 do jj=1,3 269 sscart(ii,jj)=0.0_dp 270 do kk=1,3 271 do ll=1,3 272 ! Use inverse of symop matrix here to get back to ieqqpt (inv+transpose is in symrec and in gprimd) 273 ! sscart(ii,jj)=sscart(ii,jj)+rprimd(ii,kk)*symrec(ll,kk,isym)*gprimd(ll,jj) 274 sscart(ii,jj)=sscart(ii,jj)+rprimd(ii,kk)*symrec(ll,kk,isym)*gprimd(ll,jj) 275 end do 276 end do 277 end do 278 end do 279 if (ep_scalprod==1) then 280 do ii=1,3 281 do jj=1,3 282 ss(ii,jj)=0.0_dp 283 do kk=1,3 284 do ll=1,3 285 ! Use inverse of symop matrix here to get back to ieqqpt (inv+transpose is in symrec and in gprimd) 286 ss(ii,jj)=ss(ii,jj)+rprimd(ii,kk)*symrec(ll,kk,isym)*gprimd(ll,jj) 287 end do 288 end do 289 end do 290 end do 291 else 292 do ii=1,3 293 do jj=1,3 294 ss(ii,jj) = symrel(jj,ii,isym) 295 end do 296 end do 297 end if 298 299 ss_allatoms(:,:) = zero 300 do iatom=1,crystal%natom 301 ancestor_iatom = crystal%indsym(4,isym,iatom) 302 ss_allatoms((ancestor_iatom-1)*3+1:(ancestor_iatom-1)*3+3,& 303 & (iatom-1)*3+1: (iatom-1)*3+3) = ss(1:3,1:3) 304 end do 305 306 ! ! Use inverse of symop matrix here to get back to ieqqpt 307 ! ssinv(ii,jj)=ssinv(ii,jj)+gprimd(ii,kk)*rprimd(jj,ll)*symrel(kk,ll,isym) 308 309 do isppol=1,nsppol 310 do itensor = 1, 9 311 ! multiply by the ss^{-1} matrices 312 tmp_mat2(:,:,:) = zero 313 tmp_mat(:,:,:) = reshape(gkk_qpt_new(:,itensor,:,isppol),& 314 & (/2,nbranch,nbranch/)) 315 316 317 call DGEMM ('N','N',nbranch,nbranch,nbranch,& 318 & one,ss_allatoms,nbranch,tmp_mat(1,:,:),nbranch,zero,& 319 & tmp_mat2(1,:,:),nbranch) 320 call DGEMM ('N','N',nbranch,nbranch,nbranch,& 321 & one,ss_allatoms,nbranch,tmp_mat(2,:,:),nbranch,zero,& 322 & tmp_mat2(2,:,:),nbranch) 323 324 call DGEMM ('N','T',nbranch,nbranch,nbranch,& 325 & one,tmp_mat2(1,:,:),nbranch,ss_allatoms,nbranch,zero,& 326 & tmp_mat(1,:,:),nbranch) 327 call DGEMM ('N','T',nbranch,nbranch,nbranch,& 328 & one,tmp_mat2(2,:,:),nbranch,ss_allatoms,nbranch,zero,& 329 & tmp_mat(2,:,:),nbranch) 330 331 332 gkk_qpt_tmp(:,itensor,:,isppol) = reshape (tmp_mat, (/2,nbranch*nbranch/)) 333 end do ! itensor 334 335 ! for each cartesian direction/phonon mode, rotate the tensor components 336 do imode = 1, nbranch*nbranch 337 tmp_tensor2(:,:,:) = zero 338 tmp_tensor(:,:,:) = reshape(gkk_qpt_tmp(:,:,imode,isppol),& 339 & (/2,3,3/)) 340 call DGEMM ('N','N',3,3,3,& 341 & one,sscart,3,tmp_tensor(1,:,:),3,zero,& 342 & tmp_tensor2(1,:,:),3) 343 call DGEMM ('N','T',3,3,3,& 344 & one,tmp_tensor2(1,:,:),3,sscart,3,zero,& 345 & tmp_tensor(1,:,:),3) 346 347 call DGEMM ('N','N',3,3,3,& 348 & one,sscart,3,tmp_tensor(2,:,:),3,zero,& 349 & tmp_tensor2(2,:,:),3) 350 call DGEMM ('N','T',3,3,3,& 351 & one,tmp_tensor2(2,:,:),3,sscart,3,zero,& 352 & tmp_tensor(2,:,:),3) 353 354 ! gkk_qpt_new(:,:,imode,isppol) = reshape (tmp_tensor, (/2,9/)) ! Modified by BX 355 gkk_qpt_tmp(:,:,imode,isppol) = reshape (tmp_tensor, (/2,9/)) ! Modified by BX 356 end do ! imode 357 358 if (gkk_flag (1,1,isppol,ieqqpt) == -1) then 359 gkk_flag (:,:,isppol,ieqqpt) = 0 360 end if 361 362 end do ! end isppol do 363 364 365 ! save symmetrized matrices for qpt ieqqpt 366 gamma_qpt_tr(:,:,:,:,ieqqpt) = gkk_qpt_tmp(:,:,:,:) 367 368 symmetrized_qpt(ieqqpt) = 1 369 370 end do ! end isym do 371 end do ! end itim do 372 373 end do 374 !end iqpt do 375 376 ABI_DEALLOCATE(gkk_qpt_new) 377 ABI_DEALLOCATE(gkk_qpt_tmp) 378 379 end subroutine complete_gamma_tr