TABLE OF CONTENTS


ABINIT/complete_gamma_tr [ Functions ]

[ Top ] [ 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