Line data Source code
1 : !!****m* ABINIT/m_iogkk
2 : !! NAME
3 : !! m_iogkk
4 : !!
5 : !! FUNCTION
6 : !! IO routines for GKK files
7 : !!
8 : !! COPYRIGHT
9 : !! Copyright (C) 2008-2026 ABINIT group (MVer)
10 : !! This file is distributed under the terms of the
11 : !! GNU General Public License, see ~abinit/COPYING
12 : !! or http://www.gnu.org/copyleft/gpl.txt .
13 : !!
14 : !! SOURCE
15 :
16 : #if defined HAVE_CONFIG_H
17 : #include "config.h"
18 : #endif
19 :
20 : #include "abi_common.h"
21 :
22 : module m_iogkk
23 :
24 : use defs_basis
25 : use defs_elphon
26 : use m_errors
27 : use m_abicore
28 : use m_xmpi
29 : use m_krank
30 : use m_hdr
31 :
32 : use defs_abitypes, only : MPI_type
33 : use m_numeric_tools, only : wrap2_pmhalf
34 : use m_io_tools, only : open_file, get_unit
35 : use m_matrix, only : mati3inv
36 : use m_symtk, only : littlegroup_q
37 : use m_geometry, only : phdispl_cart2red, littlegroup_pert
38 : use m_crystal, only : crystal_t
39 : use m_ifc, only : ifc_type
40 : use m_dynmat, only : d2sym3
41 : use m_ebands, only : ebands_t
42 :
43 : implicit none
44 :
45 : private
46 : !!***
47 :
48 : public :: read_gkk
49 : public :: outgkk
50 : public :: read_el_veloc
51 : !!***
52 :
53 : contains
54 : !!***
55 :
56 : !!****f* m_iogkk/read_gkk
57 : !!
58 : !! NAME
59 : !! read_gkk
60 : !!
61 : !! FUNCTION
62 : !! This routine reads in elphon matrix elements and completes them
63 : !! using the appropriate symmetries
64 : !!
65 : !! INPUTS
66 : !! elph_ds = datastructure containing elphon matrix elements
67 : !! Cryst<crystal_t>=Info on the crystal unit cell.
68 : !! Ifc<ifc_type>=Object containing the interatomic force constants.
69 : !! FSfullpqtofull = mapping of k+q to k
70 : !! n1wf = number of 1WF files to be read and analyzed
71 : !! nband = number of bands per kpoint
72 : !! unitgkk = unit of GKK file for reading
73 : !!
74 : !! OUTPUT
75 : !! elph_ds = modified gkq
76 : !! gkk_qpt = el-ph matrix elements for irreducible qpoints and
77 : !! kpoints (as a function of the reduced symmetry for the qpoint)
78 : !! gkk_flag = flag array:
79 : !! -1 -> element is missing
80 : !! 0 -> element is from symmetric qpt (Now done in complete_gkk)
81 : !! 1 -> element is from symmetric pert
82 : !! 2 -> element is kptsym of gkk file
83 : !! 3 -> element was read from gkk file
84 : !!
85 : !! SOURCE
86 :
87 15 : subroutine read_gkk(elph_ds,Cryst,ifc,Bst,FSfullpqtofull,gkk_flag,n1wf,nband,ep_prt_yambo,unitgkk)
88 :
89 : !Arguments ------------------------------------
90 : !scalars
91 : integer,intent(in) :: n1wf,nband,unitgkk,ep_prt_yambo
92 : type(crystal_t),intent(in) :: Cryst
93 : type(ifc_type),intent(in) :: ifc
94 : type(ebands_t),intent(in) :: Bst
95 : type(elph_type),intent(inout) :: elph_ds
96 : !arrays
97 : integer,intent(in) :: FSfullpqtofull(elph_ds%k_phon%nkpt,elph_ds%nqpt_full)
98 : integer,intent(out) :: gkk_flag(elph_ds%nbranch,elph_ds%nbranch,elph_ds%k_phon%my_nkpt,elph_ds%nsppol,elph_ds%nqpt_full)
99 :
100 : !Local variables-------------------------------
101 : !scalars
102 : integer :: nsppol,nbranch,nFSband,minFSband,comm,use_sym
103 : integer :: fform,i1wf,ikpt_phon,iatom1,iatom2
104 : integer :: ib,ib1,ib2,ibb,ibranch,idir,idir1,idir2,ierr,ii,ikpt1
105 : integer :: ipert,ipert1,ipert2,iqptirred,iqptfull,isppol,isym1
106 : integer :: itim1,jkpt_phon,new
107 : integer :: nsym1,qtimrev,syuse
108 : integer :: tdonecompl,test_flag,verify
109 : integer :: nqptirred_local
110 : integer :: master, me
111 : integer :: symrankkpt, ikpt1_phon, ik_this_proc
112 : real(dp) :: res,ss,timsign
113 : character(len=500) :: msg
114 15 : type(hdr_type) :: hdr1
115 : !arrays
116 30 : integer :: FSirrtok(3,elph_ds%k_phon%nkpt)
117 30 : integer :: symaf1(Cryst%nsym),symq(4,2,Cryst%nsym)
118 30 : integer :: symrc1(3,3,Cryst%nsym),symrl1(3,3,Cryst%nsym)
119 30 : integer :: tmpflg(3,Cryst%natom+2,3,Cryst%natom+2)
120 30 : real(dp) :: displ_cart(2,3*Cryst%natom,3*Cryst%natom)
121 30 : real(dp) :: displ_red(2,3*Cryst%natom,3*Cryst%natom)
122 45 : real(dp) :: eigvec(2,3*Cryst%natom,3*Cryst%natom),kpt(3),phfrq_tmp(3*Cryst%natom),redkpt(3)
123 30 : real(dp) :: qptirred_local(3,n1wf)
124 30 : real(dp) :: tnons1(3,Cryst%nsym)
125 15 : real(dp),allocatable :: eigen1(:,:,:),gkk_qpt_tmp(:,:,:,:)
126 15 : real(dp),allocatable :: h1_mat_el(:,:,:,:,:),h1_mat_el_sq(:,:,:,:,:)
127 15 : real(dp),allocatable :: qdata(:,:,:),qdata_tmp(:,:,:,:)
128 :
129 : ! *************************************************************************
130 :
131 : ABI_UNUSED(Bst%bantot)
132 :
133 15 : use_sym = 1
134 15 : nsppol = elph_ds%nsppol
135 15 : nbranch = elph_ds%nbranch
136 15 : nFSband = elph_ds%nFSband
137 15 : minFSband = elph_ds%minFSband
138 :
139 : !init values for parallelization
140 15 : comm = xmpi_world
141 15 : me = xmpi_comm_rank(comm)
142 15 : master = 0
143 :
144 90 : ABI_MALLOC_OR_DIE(h1_mat_el,(2, nFSband**2, nbranch, elph_ds%k_phon%my_nkpt, nsppol), ierr)
145 572815 : h1_mat_el= zero
146 :
147 90 : ABI_MALLOC_OR_DIE(h1_mat_el_sq,(2, nFSband**2, nbranch**2,elph_ds%k_phon%my_nkpt, nsppol), ierr)
148 2174527 : h1_mat_el_sq = zero
149 :
150 45 : ABI_MALLOC(elph_ds%qirredtofull,(elph_ds%nqptirred))
151 :
152 : !MG array to store the e-ph quantities calculated over the input Q-grid
153 90 : ABI_MALLOC(qdata_tmp,(elph_ds%nqptirred,nbranch,nsppol,3))
154 846 : qdata_tmp=zero
155 :
156 15 : nqptirred_local=0 !zero number of irred q-points found
157 735 : qptirred_local(:,:)=zero
158 :
159 172103 : gkk_flag = -1
160 :
161 15 : if (elph_ds%gkqwrite ==0) then
162 4588187 : elph_ds%gkk_qpt = zero
163 :
164 1 : else if (elph_ds%gkqwrite == 1) then
165 5 : ABI_MALLOC_OR_DIE(gkk_qpt_tmp,(2,elph_ds%ngkkband**2,nbranch**2,nsppol), ierr)
166 38 : gkk_qpt_tmp = zero
167 193 : do iqptirred=1,elph_ds%nqptirred*elph_ds%k_phon%nkpt
168 193 : write (elph_ds%unitgkq,REC=iqptirred) gkk_qpt_tmp
169 : end do
170 1 : ABI_FREE(gkk_qpt_tmp)
171 :
172 : else
173 0 : write (msg,'(a,i0)')' Wrong values for gkqwrite = ',elph_ds%gkqwrite
174 0 : ABI_BUG(msg)
175 : end if !gkqwrite
176 :
177 : !===========================================================
178 : !Loop over all files we have
179 : !read in header for perturbation
180 : !should check that all files are complete, have same header
181 : !(taking into account the symmetries for the qpoint),
182 : !represent the correct qpoints ...
183 : !MG: this task should be performed in mrggkk
184 : !===========================================================
185 :
186 60 : ABI_MALLOC(eigen1,(2,nband,nband))
187 195 : do i1wf=1,n1wf
188 :
189 180 : if (master == me) then
190 180 : write (msg,'(2a,i4,a,i4)')ch10,' read_gkk : reading 1WF header # ',i1wf,' /',n1wf
191 180 : call wrtout(std_out,msg,'COLL')
192 :
193 : ! Could check for compatibility of natom, kpt grids, ecut, qpt with DDB grid...
194 : ! MG: Also this task should be done in mrggkk
195 :
196 180 : call hdr1%fort_read(unitgkk, fform)
197 180 : if (fform == 0) then
198 0 : write (msg,'(a,i0,a)')' 1WF header number ',i1wf,' was mis-read. fform == 0'
199 0 : ABI_ERROR(msg)
200 : end if
201 :
202 180 : write(msg,'(a,i4)')' read_gkk : have read 1WF header #',i1wf
203 180 : call wrtout(std_out,msg,'COLL')
204 180 : write (msg,'(2a,i4,a)')ch10,' read_gkk : # of kpt for this perturbation: ',hdr1%nkpt,ch10
205 180 : call wrtout(std_out,msg,'COLL')
206 :
207 : end if
208 :
209 : ! broadcast data to all nodes:
210 180 : call hdr1%bcast(master, me, comm)
211 :
212 : ! Find qpoint in full grid
213 180 : new=1
214 504 : do iqptfull=1,elph_ds%nqpt_full
215 2016 : kpt(:) = hdr1%qptn(:) - elph_ds%qpt_full(:,iqptfull)
216 504 : call wrap2_pmhalf(kpt(1),redkpt(1),res)
217 504 : call wrap2_pmhalf(kpt(2),redkpt(2),res)
218 504 : call wrap2_pmhalf(kpt(3),redkpt(3),res)
219 504 : ss=redkpt(1)**2+redkpt(2)**2+redkpt(3)**2
220 504 : if(ss < tol6) then
221 : new = 0
222 : exit !exit with iqptfull
223 : end if
224 : end do !iqptfull
225 :
226 180 : if (new == 1) then
227 : ! Test should be at the end: dont care if there are additional
228 : ! qpts in gkk file which are not on the main grid. Ignore them.
229 0 : write (msg,'(4a,3es16.6,2a)')ch10,&
230 0 : & ' read_gkk : WARNING- ',ch10,&
231 0 : & ' qpoint = ',hdr1%qptn(:),ch10,&
232 0 : & ' not found in the input q-grid. Ignoring this point '
233 0 : call wrtout(ab_out,msg,'COLL')
234 0 : call wrtout(std_out,msg,'COLL')
235 0 : if (me == master) then
236 0 : do isppol=1,hdr1%nsppol
237 0 : do ikpt1=1,hdr1%nkpt
238 0 : read(unitgkk) ((eigen1(:,ii,ib),ii=1,nband),ib=1,nband)
239 : end do
240 : end do
241 : end if
242 :
243 : cycle !cycle the loop on i1wf
244 : end if !end if (new ==1)
245 :
246 :
247 : ! Check whether other pieces of the DDB have used this qpt already
248 180 : new=1
249 396 : do iqptirred=1,nqptirred_local
250 1392 : kpt(:) = qptirred_local(:,iqptirred) - hdr1%qptn(:)
251 348 : call wrap2_pmhalf(kpt(1),redkpt(1),res)
252 348 : call wrap2_pmhalf(kpt(2),redkpt(2),res)
253 348 : call wrap2_pmhalf(kpt(3),redkpt(3),res)
254 348 : ss=redkpt(1)**2+redkpt(2)**2+redkpt(3)**2
255 396 : if(ss < tol6) then
256 : new=0
257 : exit !MG We can use this information to avoid recalculating the dynamical matrix
258 : end if !but we need to use a fixed format in GKK!
259 : end do !iqptirred
260 :
261 180 : if (new==1) then !we have a new valid irreducible qpoint, add it!
262 48 : nqptirred_local = nqptirred_local+1
263 48 : if (nqptirred_local > elph_ds%nqptirred) then
264 : write (msg, '(a,a,a,i6,i6)') &
265 0 : & 'found too many qpoints in GKK file wrt anaddb input ', ch10, &
266 0 : & 'nqpt_anaddb nqpt_gkk = ', elph_ds%nqptirred, nqptirred_local
267 0 : ABI_ERROR(msg)
268 : end if
269 192 : qptirred_local(:,nqptirred_local) = hdr1%qptn(:)
270 48 : iqptirred = nqptirred_local
271 48 : tdonecompl = 0
272 1871523 : h1_mat_el = zero
273 : end if
274 :
275 : ! now iqptirred is the index of the present qpoint in the array qptirred_local
276 : ! and iqptfull is the index in the full qpt_full array for future reference
277 180 : elph_ds%qirredtofull(iqptirred) = iqptfull
278 :
279 : write (msg,'(a,i5,a,3es16.8)')&
280 180 : & ' read_gkk : full zone qpt number ',iqptfull,' is ',elph_ds%qpt_full(:,iqptfull)
281 180 : call wrtout(std_out,msg,'COLL')
282 :
283 : ! if this perturbation has already been filled (overcomplete gkk)
284 : ! check only 1st kpoint and spinpol, then check others
285 180 : verify = 0
286 180 : if (gkk_flag(hdr1%pertcase,hdr1%pertcase,1,1,elph_ds%qirredtofull(iqptirred)) /= -1) then
287 : !
288 0 : do isppol=1,nsppol
289 0 : do ik_this_proc=1,elph_ds%k_phon%my_nkpt
290 0 : if (gkk_flag(hdr1%pertcase,hdr1%pertcase,ik_this_proc,isppol,elph_ds%qirredtofull(iqptirred)) == -1) then
291 0 : write (std_out,*)" hdr1%pertcase,ik_this_proc,iqptirred",hdr1%pertcase,ik_this_proc,iqptirred
292 0 : ABI_ERROR('Partially filled perturbation ')
293 : end if
294 : end do ! ikpt_phon
295 : end do ! isppol
296 : !
297 0 : ABI_WARNING(' gkk perturbation is already filled')
298 0 : write(std_out,*)' hdr1%pertcase,iqptirred,iqptfull = ',hdr1%pertcase,iqptirred,iqptfull,&
299 0 : & gkk_flag(hdr1%pertcase,hdr1%pertcase,1,1,elph_ds%qirredtofull(iqptirred))
300 0 : verify = 1
301 0 : write (125,*) '# matrix elements for symmetric perturbation'
302 : ! Instead of reading eigen1 into void, verify == 1 checks them later on wrt values in memory
303 : end if !gkk_flag
304 :
305 : ! Examine the symmetries of the q wavevector
306 : ! these will be used to complete the perturbations for other atoms and idir
307 180 : if (ep_prt_yambo==1) then
308 : ! If one wants to print GKKs along phonon modes, it mean mixing of
309 : ! perturbations with differnt jauge. Symmetries must then be disable.
310 9 : call littlegroup_q(Cryst%nsym,qptirred_local(:,iqptirred),symq,Cryst%symrec,Cryst%symafm,qtimrev,prtvol=0,use_sym=0)
311 : else
312 171 : call littlegroup_q(Cryst%nsym,qptirred_local(:,iqptirred),symq,Cryst%symrec,Cryst%symafm,qtimrev,prtvol=0)
313 : end if
314 :
315 : ! Determine dynamical matrix, phonon frequencies and displacement vector for qpoint
316 : !call wrtout(std_out,' read_gkk: calling inpphon to calculate the dynamical matrix','COLL')
317 :
318 180 : call ifc%fourq(cryst,qptirred_local(:,iqptirred),phfrq_tmp,displ_cart,out_eigvec=eigvec)
319 :
320 : ! Get displacement vectors for all branches in reduced coordinates
321 : ! used in scalar product with H(1)_atom,idir matrix elements
322 : ! Calculate $displ_red = displ_cart \cdot gprimd$ for each phonon branch
323 :
324 180 : call phdispl_cart2red(Cryst%natom,Cryst%gprimd,displ_cart,displ_red)
325 :
326 : ! prefactors for gk+q,n\prime;k,n matrix element
327 : ! COMMENT : in decaft there is a weird term in the mass factor, of M-zval(species)
328 : ! dont know why. Not needed to reproduce decaft results, though...
329 : ! weight is squared in evaluation of
330 : ! gamma_{k,q,j} = 2 \pi omega_{q,j} sum_{nu,nu\prime} |g^{q,j}_{k+q,nu\prime; k,nu}|^2
331 : ! normally cancels with the 2 \pi omega_{q,j} factor in front of the sum...
332 :
333 : ! hdr1%pertcase = idir + (ipert-1)*3 where ipert=iatom in the interesting cases
334 180 : idir = mod (hdr1%pertcase-1,3)+1
335 180 : ipert = int(dble(hdr1%pertcase-idir)/three)+1
336 :
337 180 : write (msg,'(4a,i3,a,i3,a,i4,a)')ch10,&
338 180 : & ' read_gkk : calling littlegroup_pert to examine the symmetries of the full perturbation ',ch10,&
339 360 : & ' idir = ',idir,' ipert = ',ipert,' and Q point = ',iqptirred,ch10
340 180 : call wrtout(std_out,msg,'COLL')
341 :
342 : ! Examine the symmetries of the full perturbation these will be used to complete the kpoints
343 : ! DOESNT USE TIME REVERSAL IN littlegroup_pert except for gamma
344 :
345 180 : syuse=0
346 :
347 : call littlegroup_pert(Cryst%gprimd,idir,Cryst%indsym,ab_out,ipert,Cryst%natom,Cryst%nsym,nsym1,2,Cryst%symafm,symaf1,&
348 180 : & symq,Cryst%symrec,Cryst%symrel,symrl1,syuse,Cryst%tnons,tnons1)
349 :
350 599 : do isym1=1,nsym1
351 599 : call mati3inv(symrl1(:,:,isym1),symrc1(:,:,isym1))
352 : end do
353 56052 : FSirrtok = 0
354 :
355 : ! ========================================================
356 : ! Loop over irred kpts in file, and fill the default gkk
357 : ! ========================================================
358 :
359 : ! MG NOTE : in the present implementation, if nsppol /=1 the code stops in rchkGSheader!
360 369 : do isppol=1,hdr1%nsppol !Loop over spins is trivial? Not tested.
361 189 : write (std_out,*) ' read_gkk : isppol = ', isppol
362 :
363 14375 : do ikpt1=1,hdr1%nkpt !Loop over irred kpoints, WARNING nkpt depends on qpoint and symmetry!
364 : !
365 : ! this is the main read of the gkk matrix elements from the file (eigen1 arrays)
366 : ! it has to be done exactly nsppol*nkpt times, and the kpt_phon are completed
367 : ! where appropriate in the loop below (normally succeeding only once for each kpt)
368 : !
369 14006 : if (master == me) then
370 1551074 : read(unitgkk) ((eigen1(:,ii,ib),ii=1,nband),ib=1,nband)
371 : end if
372 :
373 : ! MPI broadcast data to all nodes:
374 14006 : call xmpi_bcast(eigen1, master, comm, ierr)
375 :
376 : ! find place of irred k in k_phon
377 : ! the kpoints in the file (kptns) could be ordered arbitrarily
378 56024 : symrankkpt = elph_ds%k_phon%krank%get_rank (hdr1%kptns(:,ikpt1)-qptirred_local(:,iqptirred))
379 14006 : ikpt1_phon = elph_ds%k_phon%krank%invrank(symrankkpt)
380 14006 : if (ikpt1_phon < 0) then
381 0 : write (msg,'(a,3es16.6,a)')' irred k ',hdr1%kptns(:,ikpt1),' was not found in full grid'
382 0 : ABI_ERROR(msg)
383 : end if
384 : ! find correspondence between this kpt_phon and the others
385 : ! symrc1 conserves perturbation as well as qpoint
386 : ! add to FSirrtok list
387 47006 : do isym1=1,nsym1
388 96350 : do itim1=0,qtimrev
389 49344 : timsign=one-two*itim1
390 1430976 : kpt(:) = timsign*matmul(symrc1(:,:,isym1), elph_ds%k_phon%kpt(:,ikpt1_phon))
391 :
392 49344 : symrankkpt = elph_ds%k_phon%krank%get_rank (kpt)
393 49344 : jkpt_phon = elph_ds%k_phon%krank%invrank(symrankkpt)
394 :
395 82344 : if (jkpt_phon > 0) then
396 49344 : FSirrtok(1,jkpt_phon) = ikpt1_phon
397 49344 : FSirrtok(2,jkpt_phon) = isym1
398 49344 : FSirrtok(3,jkpt_phon) = itim1
399 : else
400 : write (msg,'(a,3es16.6,a,i5,a,i4,a)')&
401 0 : & ' sym equivalent of kpt ',hdr1%kptns(:,ikpt1),' by sym ',&
402 0 : & isym1,' and itime ',itim1,' was not found'
403 0 : ABI_ERROR(msg)
404 : end if
405 : end do !itim1
406 : end do !isim1
407 :
408 :
409 : !
410 : ! Here check if the symmetry-copied gkk at new k point is equal to the one found in the file for non-irreducible point
411 : ! NB This is DEBUG code
412 : !
413 14006 : if (verify == 1 .and. elph_ds%k_phon%my_kpt(ikpt1_phon) == me) then
414 0 : do ik_this_proc = 1, elph_ds%k_phon%my_nkpt
415 0 : if (elph_ds%k_phon%my_ikpt(ik_this_proc) == ikpt1_phon) exit
416 : end do
417 0 : do ib1=1,nFSband
418 0 : do ib2=1,nFSband
419 0 : ibb = (ib1-1)*nFSband+ib2
420 0 : write (125,'(2(2E16.6,2x))') h1_mat_el(:,ibb,hdr1%pertcase,ik_this_proc,isppol),&
421 0 : & eigen1(:,minFSband-1+ib2,minFSband-1+ib1)
422 : end do
423 : end do
424 : end if !verify end DEBUG code
425 :
426 :
427 1786457 : do ik_this_proc = 1, elph_ds%k_phon%my_nkpt
428 : ! should I be dealing with this k-point?
429 1758256 : jkpt_phon = elph_ds%k_phon%my_ikpt(ik_this_proc)
430 :
431 : ! does present ikpt1 contribute to this k-point?
432 1758256 : if (FSirrtok(1,jkpt_phon) /= ikpt1_phon) cycle
433 :
434 : ! if this kpoint has already been filled (overcomplete gkk)
435 36428 : if (gkk_flag(hdr1%pertcase,hdr1%pertcase,ik_this_proc,isppol,elph_ds%qirredtofull(iqptirred)) /= -1) then
436 22388 : ABI_WARNING("gkk element is already filled")
437 22388 : write(std_out,*)' hdr1%pertcase,ik_this_proc,isppol,iqptirred = ',&
438 22388 : & hdr1%pertcase,ik_this_proc,isppol,iqptirred,&
439 44776 : & gkk_flag(hdr1%pertcase,hdr1%pertcase,ik_this_proc,isppol,elph_ds%qirredtofull(iqptirred))
440 : ! exit
441 : end if !gkk_flag
442 :
443 : ! ===============================================================
444 : ! TODO: if there is a phase factor in swapping k-points, insert it here in copy to h1_mat_el
445 : ! as a function of symops in FSirrtok
446 : ! complete gkk for symmetric ikpt_phon with sym1 which conserve
447 : ! the full perturbation+qpoint
448 : ! Not tested explicitly, but the results for Pb using reduced kpts look good
449 : ! should do same RF calculation with nsym=1 and check
450 : ! ===============================================================
451 :
452 : ! save this kpoint
453 269924 : do ib1=1,nFSband
454 1802356 : do ib2=1,nFSband
455 1532432 : ibb = (ib1-1)*nFSband+ib2
456 :
457 : ! real
458 1532432 : res=eigen1(1,minFSband-1+ib2,minFSband-1+ib1)
459 1532432 : h1_mat_el(1,ibb,hdr1%pertcase,ik_this_proc,isppol) = res
460 :
461 : ! imag
462 1532432 : res=eigen1(2,minFSband-1+ib2,minFSband-1+ib1)
463 1765928 : h1_mat_el(2,ibb,hdr1%pertcase,ik_this_proc,isppol) = res
464 : end do !ib2
465 : end do !ib1
466 : ! if jkpt is equal to ikpt1_phon (if clause above) flag == 3
467 50434 : if (FSirrtok(2,jkpt_phon) == 1) then
468 11944 : gkk_flag(hdr1%pertcase,hdr1%pertcase,ik_this_proc,isppol,elph_ds%qirredtofull(iqptirred)) = 3
469 : ! if jkpt_phon comes from ikpt1_phon flag == 2 with some symop
470 : else
471 24484 : gkk_flag(hdr1%pertcase,hdr1%pertcase,ik_this_proc,isppol,elph_ds%qirredtofull(iqptirred)) = 2
472 : end if
473 :
474 : end do !jkpt_phon
475 :
476 : ! ===============================================================
477 : ! we now have contribution to g(k+q,k; \kappa,\alpha) from one
478 : ! kpoint,and one perturbation,
479 : ! NB: each perturbation will contribute to all the modes later!
480 : !
481 : ! SHOULD ONLY DO THIS FOR THE SYMS NEEDED
482 : ! TO COMPLETE THE PERTURBATIONS!!!
483 : ! ================================================================
484 :
485 : end do !ikpt1
486 : end do !isppol
487 :
488 : ! 14 Jan 2014 removed test on verify - in new scheme full BZ is read in and should be used to avoid phase errors
489 : ! if (verify == 1) cycle
490 :
491 : ! Checks on irred grid provided and on gkk_flag accumulated up to now
492 180 : if (elph_ds%tuniformgrid == 1) then ! check if irred kpoints found reconstitute the FS kpts
493 14148 : do ikpt_phon=1,elph_ds%k_phon%nkpt
494 14148 : if (FSirrtok(1,ikpt_phon) == 0) then
495 : write(msg,'(a,3es16.6,2a)')&
496 0 : & ' kpt = ',elph_ds%k_phon%kpt(:,ikpt_phon),ch10,&
497 0 : & ' is not the symmetric of one of those found in the GKK file'
498 0 : ABI_ERROR(msg)
499 : end if
500 : end do !ikpt_phon
501 :
502 : ! normally at this point we have used all the gkk for all kpoints on the FS
503 : ! for the given irred perturbation: check
504 14148 : do ik_this_proc = 1, elph_ds%k_phon%my_nkpt
505 13968 : ikpt_phon = elph_ds%k_phon%my_ikpt(ik_this_proc)
506 :
507 14148 : if (gkk_flag(hdr1%pertcase, hdr1%pertcase, ik_this_proc, 1, elph_ds%qirredtofull(iqptirred)) == -1) then
508 : write (msg,'(a,i3,a,3es18.6,2a,i3,a,i3,a,3es18.6,a,a,i4,a,a)')&
509 0 : & ' For irreducible qpt ', iqptirred,' = ',qptirred_local(:,iqptirred),ch10, &
510 0 : & ' the gkk element : pertcase = ',hdr1%pertcase,' ik_this_proc = ',ik_this_proc, &
511 0 : & ' kpt = ',elph_ds%k_phon%kpt(:,ikpt_phon),ch10,&
512 0 : & ' and isppol ',1,ch10,&
513 0 : & ' was not found by symmetry operations on the irreducible kpoints given'
514 0 : ABI_ERROR(msg)
515 : end if
516 : end do !ikpt_phon
517 : end if ! end elph_ds%tuniformgrid == 1 checks
518 :
519 180 : write(msg,'(a,i0)')' read_gkk : Done completing the kpoints for pertcase ',hdr1%pertcase
520 180 : call wrtout(std_out,msg,'COLL')
521 :
522 28116 : tmpflg(:,:,:,:) = 0
523 :
524 720 : do idir1=1,3
525 1476 : do iatom1=1,Cryst%natom
526 756 : ipert1 = (iatom1-1)*3+idir1
527 3564 : do idir2=1,3
528 6588 : do iatom2=1,Cryst%natom
529 3564 : ipert2 = (iatom2-1)*3+idir2
530 3564 : if (gkk_flag(ipert1,ipert1,1,1,elph_ds%qirredtofull(iqptirred)) >= 0 .and. &
531 2268 : & gkk_flag(ipert2,ipert2,1,1,elph_ds%qirredtofull(iqptirred)) >= 0) then
532 1596 : tmpflg(idir1,iatom1,idir2,iatom2) = 1
533 : end if
534 : end do
535 : end do
536 : end do
537 : end do
538 :
539 :
540 : ! ===============================================
541 : ! Full test: need all perturbations explicitly
542 : ! ===============================================
543 :
544 180 : test_flag = 0
545 5940 : if (sum(tmpflg(:,1:Cryst%natom,:,1:Cryst%natom)) == (3*Cryst%natom)**2 .and. tdonecompl == 0) test_flag = 1
546 :
547 180 : write(std_out,*)'read_gkk: tdonecompl = ', tdonecompl
548 :
549 : ! de-activate completion of perts by symmetry for now.
550 : ! Must be called when all irreducible perturbations are in memory!!!!
551 180 : if (test_flag == 1 .and. tdonecompl == 0) then
552 :
553 : ! write(std_out,*) ' read_gkk : enter fxgkkphase before completeperts'
554 : ! call fxgkkphase(elph_ds,gkk_flag,h1_mat_el,iqptirred)
555 :
556 48 : if (ep_prt_yambo==1) then
557 3 : if (elph_ds%k_phon%my_nkpt /= elph_ds%k_phon%nkpt) then
558 0 : write (msg, '(a)') 'prt_gkk_yambo can not handle parallel anaddb yet'
559 0 : ABI_ERROR(msg)
560 : end if
561 : call prt_gkk_yambo(displ_cart,displ_red,elph_ds%k_phon%kpt,h1_mat_el,iqptirred,&
562 3 : & Cryst%natom,nFSband,elph_ds%k_phon%my_nkpt,phfrq_tmp,hdr1%qptn)
563 : end if
564 :
565 : ! ========================================================================
566 : ! Now use more general symops to complete the other equivalent
567 : ! perturbations: the kpoints are also shuffled by these symops
568 : ! afterwards h1_mat_el_sq contains gamma_\tau\alpha,\tau'\alpha' in reduced coordinates
569 : !
570 : ! \gamma_{\tau'\alpha',\tau\alpha} =
571 : ! <psi_{k+q,ib2}| H(1)_{\tau'\alpha'}| psi_{k,ib1}>* \cdot
572 : ! <psi_{k+q,ib2}| H(1)_{\tau \alpha }| psi_{k,ib1}>
573 : !
574 : ! ========================================================================
575 :
576 : call completeperts(Cryst,nbranch,nFSband,elph_ds%k_phon%my_nkpt,nsppol,&
577 48 : & gkk_flag(:,:,:,:,elph_ds%qirredtofull(iqptirred)),h1_mat_el,h1_mat_el_sq,qptirred_local(:,iqptirred),symq,qtimrev)
578 :
579 48 : tdonecompl = 1
580 : end if
581 :
582 : ! ==============================================================
583 : ! if we have all the perturbations for this qpoint, proceed
584 : ! with scalar product, norm squared, and add weight factors
585 : !
586 : ! SHOULD HAVE A TEST SO h1_mat_el IS NOT OVERWRITTEN
587 : ! BEFORE PREVIOUS QPOINT IS FINISHED!!!!!
588 : ! ==============================================================
589 :
590 180 : test_flag = 1
591 369 : do isppol=1,nsppol
592 14409 : do ik_this_proc = 1, elph_ds%k_phon%my_nkpt
593 47493 : do ibranch=1,nbranch
594 47304 : if (gkk_flag (ibranch,ibranch,ik_this_proc,isppol,elph_ds%qirredtofull(iqptirred)) == -1) then
595 : test_flag = 0
596 : exit
597 : end if
598 : end do
599 : end do
600 : end do
601 :
602 180 : if (test_flag /= 0) then
603 48 : call wrtout(std_out,' read_gkk : enter normsq_gkq',"COLL")
604 :
605 : ! MG temporary array to save ph-linewidths before Fourier interpolation
606 240 : ABI_MALLOC(qdata,(nbranch,nsppol,3))
607 912 : qdata(:,:,:)=zero
608 :
609 : call normsq_gkq(displ_red,eigvec,elph_ds,FSfullpqtofull,&
610 48 : & h1_mat_el_sq,iqptirred,phfrq_tmp,qptirred_local,qdata)
611 :
612 : ! save gkk_qpt, eventually to disk, for bands up to ngkkband,
613 : ! NB: if the sum over bands has been performed ngkkband is 1 instead of nFSband
614 48 : if (elph_ds%gkqwrite == 0) then
615 4588173 : elph_ds%gkk_qpt(:,:,:,:,:,iqptirred) = h1_mat_el_sq(:,1:elph_ds%ngkkband*elph_ds%ngkkband,:,:,:)
616 : else
617 : ! write all kpoints to disk
618 3 : write (std_out,*) 'size of record to be written: ', 8 * 2*elph_ds%ngkkband*elph_ds%ngkkband*&
619 6 : & elph_ds%nbranch*elph_ds%nbranch*elph_ds%k_phon%my_nkpt*elph_ds%nsppol
620 3 : inquire(unit=elph_ds%unitgkq, recl=isppol)
621 3 : write (std_out,*) 'recl =', isppol
622 3 : write (std_out,*) 'iqptirred ', iqptirred
623 195 : do ik_this_proc = 1, elph_ds%k_phon%my_nkpt
624 : write (elph_ds%unitgkq,REC=((iqptirred-1)*elph_ds%k_phon%my_nkpt+ik_this_proc)) &
625 195 : & h1_mat_el_sq(:,1:elph_ds%ngkkband*elph_ds%ngkkband,:,ik_this_proc,:)
626 : end do
627 : end if
628 :
629 912 : qdata_tmp(iqptirred,:,:,:)=qdata(:,:,:)
630 48 : ABI_FREE(qdata)
631 : end if
632 :
633 195 : call hdr1%free()
634 :
635 : end do !of i1wf
636 :
637 : !got all the gkk perturbations
638 :
639 15 : ABI_FREE(eigen1)
640 15 : ABI_FREE(h1_mat_el)
641 15 : ABI_FREE(h1_mat_el_sq)
642 :
643 15 : if (nqptirred_local /= elph_ds%nqptirred) then
644 : write (msg, '(3a,i0,i0)') &
645 0 : & ' Found wrong number of qpoints in GKK file wrt anaddb input ', ch10, &
646 0 : & ' nqpt_anaddb nqpt_gkk = ', elph_ds%nqptirred, nqptirred_local
647 0 : ABI_ERROR(msg)
648 : end if
649 :
650 : !normally at this point we have the gkk for all kpoints on the FS
651 : !for all the perturbations. Otherwise a 1WF file is missing.
652 : !NOTE: still havent checked the qpoint grid completeness
653 63 : do iqptirred=1,elph_ds%nqptirred
654 114 : do isppol=1,nsppol
655 4203 : do ik_this_proc = 1, elph_ds%k_phon%my_nkpt
656 4104 : ikpt_phon = elph_ds%k_phon%my_ikpt(ik_this_proc)
657 18195 : do ipert=1,nbranch
658 18144 : if (gkk_flag(ipert,ipert,ik_this_proc,isppol,elph_ds%qirredtofull(iqptirred)) == -1) then
659 : write (msg,'(a,i5,1x,i5,1x,i5,1x,i5,a,a)')&
660 0 : & ' gkk element',ipert,ikpt_phon,isppol,iqptirred,' was not found by symmetry operations ',&
661 0 : & ' on the irreducible perturbations and qpoints given'
662 0 : ABI_ERROR(msg)
663 : end if
664 : end do !ipert
665 : end do !ik_this_proc
666 : end do !isppol
667 : end do !iqptirred
668 :
669 15 : call wrtout(std_out,'read_gkk : done completing the perturbations (and checked!)','COLL')
670 :
671 : !MG save phonon frequencies, ph-linewidths and lambda(q,n) values before Fourier interpolation
672 90 : ABI_MALLOC(elph_ds%qgrid_data,(elph_ds%nqptirred,nbranch,nsppol,3))
673 :
674 63 : do iqptirred=1,elph_ds%nqptirred
675 927 : elph_ds%qgrid_data(iqptirred,:,:,:)=qdata_tmp(iqptirred,:,:,:)
676 : end do
677 :
678 15 : ABI_FREE(qdata_tmp)
679 :
680 15 : end subroutine read_gkk
681 : !!***
682 :
683 : !!****f* m_iogkk/outgkk
684 : !! NAME
685 : !! outgkk
686 : !!
687 : !! FUNCTION
688 : !! output gkk file for one perturbation (used for elphon calculations in anaddb)
689 : !!
690 : !! INPUTS
691 : !! bantot0 = total number of bands for all kpoints
692 : !! bantot1 = total number of matrix elements for 1st order eigenvalues
693 : !! eigen0 = GS eigenvalues
694 : !! eigen1 = response function 1st order eigenvalue matrix
695 : !! hdr0 = GS header
696 : !! hdr1 = RF header
697 : !! mpi_enreg=information about MPI parallelization
698 : !!
699 : !! SOURCE
700 :
701 90 : subroutine outgkk(bantot0,bantot1,outfile,eigen0,eigen1,hdr0,hdr1,mpi_enreg,phasecg)
702 :
703 : !Arguments ------------------------------------
704 : !scalars
705 : integer,intent(in) :: bantot0,bantot1
706 : character(len=fnlen),intent(in) :: outfile
707 : type(MPI_type),intent(in) :: mpi_enreg
708 : type(hdr_type),intent(inout) :: hdr0,hdr1
709 : !arrays
710 : real(dp),intent(in) :: eigen0(bantot0),eigen1(2*bantot1)
711 : real(dp),intent(in) :: phasecg(2,bantot1)
712 :
713 : !Local variables-------------------------------
714 : !scalars
715 : integer :: fform,iband,ikpt,isppol,me,ntot,unitout
716 : integer :: iband_off, mband, ierr
717 : character(len=500) :: msg
718 90 : real(dp), allocatable :: tmpeig(:)
719 :
720 : ! *************************************************************************
721 :
722 : !only master should be writing to disk
723 : !Init me
724 90 : me=mpi_enreg%me_kpt
725 90 : if (me /= 0) return
726 :
727 90 : call wrtout(std_out,' writing gkk file: '//outfile,"COLL")
728 :
729 : !initializations
730 90 : fform = 42
731 90 : ntot = 1
732 :
733 : !open gkk file
734 90 : if (open_file(outfile, msg, newunit=unitout, form='unformatted', status='unknown', action="write") /= 0) then
735 0 : ABI_ERROR(msg)
736 : end if
737 :
738 : !output GS header
739 90 : call hdr0%fort_write(unitout, fform, ierr)
740 90 : ABI_CHECK(ierr == 0 , "hdr_fort_write returned ierr != 0")
741 :
742 : !output GS eigenvalues
743 90 : iband=0
744 180 : do isppol=1,hdr0%nsppol
745 6468 : do ikpt=1,hdr0%nkpt
746 6288 : write (unitout) eigen0(iband+1:iband+hdr0%nband(ikpt))
747 6378 : iband=iband+hdr0%nband(ikpt)
748 : end do
749 : end do
750 :
751 : !output number of gkk in this file (1)
752 90 : write (unitout) ntot
753 :
754 : !output RF header
755 90 : call hdr1%fort_write(unitout, fform, ierr)
756 90 : ABI_CHECK(ierr == 0 , "hdr_fort_write returned ierr != 0")
757 :
758 : !output RF eigenvalues
759 6378 : mband = maxval(hdr1%nband(:))
760 270 : ABI_MALLOC(tmpeig,(2*mband**2))
761 90 : iband_off = 0
762 90 : tmpeig(1) = phasecg(1, 1)
763 180 : do isppol = 1, hdr1%nsppol
764 6468 : do ikpt = 1, hdr1%nkpt
765 1177488 : tmpeig = zero
766 591888 : do iband = 1, hdr1%nband(ikpt)**2
767 585600 : tmpeig (2*(iband-1)+1) = eigen1(2*(iband_off+iband-1)+1)
768 591888 : tmpeig (2*(iband-1)+2) = eigen1(2*(iband_off+iband-1)+2)
769 : end do
770 6288 : write (unitout) tmpeig(1:2*hdr1%nband(ikpt)**2)
771 6378 : iband_off = iband_off + hdr1%nband(ikpt)**2
772 : end do
773 : end do
774 90 : ABI_FREE(tmpeig)
775 :
776 : !close gkk file
777 90 : close (unitout)
778 :
779 : end subroutine outgkk
780 : !!***
781 :
782 : !!****f* m_iogkk/prt_gkk_yambo
783 : !!
784 : !! NAME
785 : !! prt_gkk_yambo
786 : !!
787 : !! FUNCTION
788 : !! This routine outputs el-phon related quantities for the yambo code at 1
789 : !! q-point
790 : !!
791 : !! INPUTS
792 : !! displ_cart = phonon displacement vectors for this q-point in Cartesian coordinates.
793 : !! displ_red = phonon displacement vectors for this q-point, in reduced coordinates
794 : !! elph_ds = datastructure containing elphon matrix elements
795 : !! h1_mat_el = matrix elements of first order hamiltonian for present q-point,
796 : !! all perturbations
797 : !! iqptfull = index of present q-point in full array of q-points
798 : !! irredpert = index of irreducible perturbation (atom displaced)
799 : !! natom = number of atoms
800 : !! phfrq = phonon frequencies at present q-point
801 : !! qptn = q-point we will print for
802 : !!
803 : !! OUTPUT
804 : !! only writes to a file
805 : !!
806 : !! NOTES
807 : !!
808 : !! SOURCE
809 :
810 3 : subroutine prt_gkk_yambo(displ_cart,displ_red,kpt_phon,h1_mat_el,iqpt,&
811 3 : & natom,nFSband,nkpt_phon,phfrq,qptn)
812 :
813 : !Arguments ------------------------------------
814 : !scalars
815 : integer,intent(in) :: natom,iqpt
816 : integer,intent(in) :: nFSband,nkpt_phon
817 : !arrays
818 : real(dp),intent(in) :: kpt_phon(3,nkpt_phon)
819 : real(dp),intent(in) :: h1_mat_el(2,nFSband*nFSband,3*natom,nkpt_phon,1)
820 : real(dp),intent(in) :: phfrq(3*natom)
821 : real(dp),intent(in) :: displ_cart(2,3*natom,3*natom)
822 : real(dp),intent(in) :: displ_red(2,3*natom,3*natom)
823 : real(dp),intent(in) :: qptn(3)
824 :
825 : !Local variables-------------------------------
826 : !scalars
827 : integer, save :: firsttime=1
828 : integer :: outunit,ikpt,imode,iband,ibandp,iatom,idir,ibandindex
829 : integer :: jmode, outunit2, outunit3
830 : !arrays
831 : real(dp) :: gkk_mode_dep(2)
832 : ! *************************************************************************
833 :
834 : !if first time round:
835 3 : if (firsttime==1) then
836 1 : firsttime=0
837 : ! squash file
838 1 : outunit=get_unit()
839 1 : open (unit=outunit,file="yambo_elphon_data",status="REPLACE")
840 1 : outunit2=get_unit()
841 1 : open (unit=outunit2,file="yambo_elphon_gkk_bymode",status="replace")
842 1 : outunit3=get_unit()
843 1 : open (unit=outunit3,file="yambo_elphon_gkksqtw_bymode",status="replace")
844 :
845 : ! write dimensions
846 1 : write (outunit,'(a,I6)') 'number of el atoms ', natom
847 1 : write (outunit2,'(a,I6)') 'number of el atoms ', natom
848 1 : write (outunit3,'(a,I6)') 'number of el atoms ', natom
849 1 : write (outunit,'(a,I6)') 'number of ph modes ', 3*natom
850 1 : write (outunit2,'(a,I6)') 'number of ph modes ', 3*natom
851 1 : write (outunit3,'(a,I6)') 'number of ph modes ', 3*natom
852 1 : write (outunit,'(a,I6)') 'number of el bands ', nFSband
853 1 : write (outunit2,'(a,I6)') 'number of el bands ', nFSband
854 1 : write (outunit3,'(a,I6)') 'number of el bands ', nFSband
855 :
856 : ! write k-points
857 1 : write (outunit,'(a,I6)') 'number of k-points ', nkpt_phon
858 1 : write (outunit2,'(a,I6)') 'number of k-points ', nkpt_phon
859 1 : write (outunit3,'(a,I6)') 'number of k-points ', nkpt_phon
860 65 : do ikpt=1,nkpt_phon
861 64 : write (outunit,'(a,I6,3E20.10)') 'reduced coord kpoint no ', ikpt, kpt_phon(:,ikpt)
862 64 : write (outunit2,'(a,I6,3E20.10)') 'reduced coord kpoint no ', ikpt, kpt_phon(:,ikpt)
863 65 : write (outunit3,'(a,I6,3E20.10)') 'reduced coord kpoint no ', ikpt, kpt_phon(:,ikpt)
864 : end do
865 :
866 : ! band energies are not accessible this deep in the code: simpler to get them
867 : ! from elsewhere
868 :
869 1 : close (outunit)
870 1 : close (outunit2)
871 1 : close (outunit3)
872 : end if ! first time round
873 :
874 : !open file
875 3 : outunit=get_unit()
876 3 : open (unit=outunit,file="yambo_elphon_data",status="unknown",position="append")
877 :
878 : !qpoint
879 3 : write (outunit,'(a,I6,3E20.10)') 'reduced coord qpoint no ', iqpt, qptn(:)
880 :
881 : !frequencies
882 12 : do imode=1,3*natom
883 12 : write (outunit,'(a,I6,3E20.10)') 'phonon freq no ', imode, phfrq(imode)
884 : end do
885 :
886 : !displacement vector
887 12 : do imode=1,3*natom
888 9 : write (outunit,'(a,I6,3E20.10)') 'phonon displ vec no ', imode
889 21 : do iatom=1,natom
890 18 : write (outunit,'(3(2E20.10,2x))') displ_cart(:,(iatom-1)*3+1:iatom*3,imode)
891 : end do
892 : end do
893 :
894 : !the beef: matrix elements of the first order hamiltonian for displacement of
895 : !all atoms along all reduced directions
896 3 : write (outunit,'(a)') ' matrix elements of all perturbations for this q-point'
897 195 : do ikpt=1,nkpt_phon
898 192 : write (outunit,'(a,I6)') ' kpoint ', ikpt
899 192 : imode=0
900 387 : do iatom=1,natom
901 960 : do idir=1,3
902 576 : imode=imode+1
903 576 : write (outunit,'(a,I6,I6)') ' atom, direction = ', iatom,idir
904 576 : ibandindex=0
905 6528 : do iband=1,nFSband
906 63936 : do ibandp=1,nFSband
907 57600 : ibandindex=ibandindex+1
908 57600 : write (outunit,'(a,I6,I6,2E20.10)') ' mat el for n,np ', iband,ibandp,&
909 120960 : & h1_mat_el(:,ibandindex,imode,ikpt,1)
910 : end do !bandp
911 : end do !band
912 : end do !dir
913 : end do !atom
914 : end do
915 :
916 : !blank line
917 3 : write (outunit,*)
918 3 : close (outunit)
919 :
920 3 : outunit2=get_unit()
921 3 : open (unit=outunit2,file="yambo_elphon_gkk_bymode",status="unknown",position="append")
922 3 : outunit3=get_unit()
923 3 : open (unit=outunit3,file="yambo_elphon_gkksqtw_bymode",status="unknown",position="append")
924 :
925 : !qpoint
926 3 : write (outunit2,'(a,I6,3E20.10)') 'reduced coord qpoint no ', iqpt, qptn(:)
927 3 : write (outunit3,'(a,I6,3E20.10)') 'reduced coord qpoint no ', iqpt, qptn(:)
928 :
929 : !print out mode-dependent matrix elements
930 3 : write (outunit2,'(a)') ' matrix elements of all phonon modes for this q-point'
931 3 : write (outunit3,'(a)') ' 1/w**1/2 times matrix elements of all phonon modes for this q-point'
932 195 : do ikpt=1,nkpt_phon
933 192 : write (outunit2,'(a,I6)') ' kpoint ', ikpt
934 192 : write (outunit3,'(a,I6)') ' kpoint ', ikpt
935 192 : ibandindex=0
936 2115 : do iband=1,nFSband
937 21312 : do ibandp=1,nFSband
938 19200 : ibandindex=ibandindex+1
939 19200 : write (outunit2,'(a,I6,I6)') ' el bands n,np ', iband,ibandp
940 19200 : write (outunit3,'(a,I6,I6)') ' el bands n,np ', iband,ibandp
941 78720 : do imode=1,3*natom
942 : ! gkk_mode_dep = cg_zdotc(3*natom,displ_red(:,:,imode),h1_mat_el(:,ibandindex,:,ikpt,1))
943 57600 : gkk_mode_dep = zero
944 230400 : do jmode=1,3*natom
945 : gkk_mode_dep(1) = gkk_mode_dep(1) &
946 : & + displ_red(1,jmode,imode)*h1_mat_el(1,ibandindex,jmode,ikpt,1) &
947 172800 : & + displ_red(2,jmode,imode)*h1_mat_el(2,ibandindex,jmode,ikpt,1)
948 : gkk_mode_dep(2) = gkk_mode_dep(2) &
949 : & + displ_red(1,jmode,imode)*h1_mat_el(2,ibandindex,jmode,ikpt,1) &
950 230400 : & - displ_red(2,jmode,imode)*h1_mat_el(1,ibandindex,jmode,ikpt,1)
951 : end do
952 57600 : write (outunit2,'(a,I6,2E20.10)') ' mat el for phonon mode num = ', imode, gkk_mode_dep
953 57600 : write (outunit3,'(a,I6,2E20.10)') ' 1/w**1/2 * mat el for phonon mode num = ', &
954 249600 : & imode, gkk_mode_dep/sqrt(two*abs(phfrq(imode))+tol10)
955 : end do !imode
956 : end do !bandp
957 : end do !band
958 : end do
959 : !blank line
960 3 : write (outunit2,*)
961 3 : write (outunit3,*)
962 :
963 3 : close (outunit2)
964 3 : close (outunit3)
965 :
966 3 : end subroutine prt_gkk_yambo
967 : !!***
968 :
969 : !!****f* m_iogkk/read_el_veloc
970 : !!
971 : !! NAME
972 : !! read_el_veloc
973 : !!
974 : !! FUNCTION
975 : !! This routine reads the velocities of the electronic GS
976 : !! for all kpts and bands
977 : !! then maps them into the FS kpt states
978 : !!
979 : !! COPYRIGHT
980 : !! Copyright (C) 2002-2026 ABINIT group (JPCroc) based on conducti
981 : !! This file is distributed under the terms of the
982 : !! GNU General Public License, see ~abinit/COPYING
983 : !! or http://www.gnu.org/copyleft/gpl.txt .
984 : !! For the initials of contributors, see ~abinit/doc/developers/contributors.txt .
985 : !!
986 : !! INPUTS
987 : !! nkpt_in = number of kpoints according to parent routine
988 : !! nband_in = number of bands according to parent routine
989 : !! nsppol_in = number of spin polarizations
990 : !!
991 : !! OUTPUT
992 : !! el_veloc(nkpt_in,nband_in,3)
993 : !!
994 : !! SOURCE
995 :
996 5 : subroutine read_el_veloc(nband_in,nkpt_in,kpt_in,nsppol_in,elph_tr_ds)
997 :
998 : !Arguments -----------------------------------
999 : !scalars
1000 : integer, intent(in) :: nband_in,nkpt_in,nsppol_in
1001 : type(elph_tr_type), intent(inout) :: elph_tr_ds
1002 : real(dp), intent(in) :: kpt_in(3,nkpt_in)
1003 :
1004 : !Local variables-------------------------------
1005 : !scalars
1006 : integer :: bd2tot_index, iband,ii,ikpt, ikpt_ddk, isppol,l1,mband
1007 : integer :: bantot1, unit_ddk, symrankkpt
1008 : character(len=fnlen) :: filnam1,filnam2,filnam3
1009 : character(len=500) :: msg
1010 5 : type(hdr_type) :: hdr1
1011 5 : type(krank_t) :: krank
1012 : !arrays
1013 : real(dp) :: im_el_veloc(3)
1014 5 : real(dp),allocatable :: eig1_k(:,:), eigen11(:),eigen12(:),eigen13(:)
1015 : ! *********************************************************************************
1016 :
1017 : !Read data file name
1018 : !TODO: this should be standardized and read in anaddb always, not
1019 : !conditionally. Otherwise when new files are added to the anaddb files
1020 : !file... Catastrophe!
1021 :
1022 5 : write(std_out,*)'enter read_el_veloc '
1023 :
1024 : !Read data file
1025 5 : if (open_file(elph_tr_ds%ddkfilename,msg,newunit=unit_ddk,form='formatted') /= 0) then
1026 0 : ABI_ERROR(msg)
1027 : end if
1028 :
1029 5 : rewind(unit_ddk)
1030 5 : read(unit_ddk,'(a)')filnam1 ! first ddk file
1031 5 : read(unit_ddk,'(a)')filnam2 ! second ddk file
1032 5 : read(unit_ddk,'(a)')filnam3 ! third ddk file
1033 5 : close (unit_ddk)
1034 :
1035 5 : bantot1 = 2*nband_in**2*nkpt_in*nsppol_in
1036 :
1037 5 : call inpgkk(eigen11,filnam1,hdr1)
1038 5 : call hdr1%free()
1039 :
1040 5 : call inpgkk(eigen12,filnam2,hdr1)
1041 5 : call hdr1%free()
1042 :
1043 : !we use the hdr1 from the last call - should add some consistency
1044 : !testing here, we are trusting users not to mix different ddk files...
1045 5 : call inpgkk(eigen13,filnam3,hdr1)
1046 :
1047 : !Extract info from the header
1048 5 : if(hdr1%nsppol /= nsppol_in) then
1049 0 : ABI_ERROR('nsspol /= input nsppol')
1050 : end if
1051 :
1052 : !Get mband, as the maximum value of nband(nkpt)
1053 429 : mband=maxval(hdr1%nband(1:hdr1%nkpt))
1054 5 : if (mband /= nband_in) then
1055 0 : ABI_ERROR('nband_in input to read_el_veloc is inconsistent with mband')
1056 : end if
1057 :
1058 5 : write(std_out,*)
1059 5 : write(std_out,*) 'readings from read_el_veloc header'
1060 5 : write(std_out,'(a,i8)') ' natom =',hdr1%natom
1061 5 : write(std_out,'(a,3i8)') ' nkpt,nband_in,mband =',hdr1%nkpt,nband_in,mband
1062 5 : write(std_out,'(a, f10.5,a)' ) ' ecut =',hdr1%ecut,' Ha'
1063 5 : write(std_out,'(a,e15.5,a,e15.5,a)' )' fermie =',hdr1%fermie,' Ha ',hdr1%fermie*Ha_eV,' eV'
1064 :
1065 20 : ABI_MALLOC(eig1_k,(2*nband_in**2,3))
1066 5 : bd2tot_index = 0
1067 12793 : elph_tr_ds%el_veloc=zero
1068 :
1069 : !need correspondence between the DDK kpoints and the kpt_phon
1070 5 : call krank%init(hdr1%nkpt, hdr1%kptns)
1071 :
1072 10 : do isppol=1,nsppol_in
1073 : im_el_veloc(:)=zero
1074 434 : do ikpt=1,nkpt_in
1075 424 : symrankkpt = krank%get_rank (kpt_in(:,ikpt))
1076 424 : ikpt_ddk = krank%invrank(symrankkpt)
1077 424 : if (ikpt_ddk == -1) then
1078 0 : write(std_out,*)'read_el_veloc ******** error in correspondence between ddk and gkk kpoint sets'
1079 0 : write(std_out,*)' kpt sets in gkk and ddk files must agree.'
1080 0 : ABI_ERROR("Aborting now")
1081 : end if
1082 424 : bd2tot_index=2*nband_in**2*(ikpt_ddk-1)
1083 :
1084 : ! first derivative eigenvalues for k-point
1085 84072 : eig1_k(:,1)=eigen11(1+bd2tot_index:2*nband_in**2+bd2tot_index)
1086 84072 : eig1_k(:,2)=eigen12(1+bd2tot_index:2*nband_in**2+bd2tot_index)
1087 84072 : eig1_k(:,3)=eigen13(1+bd2tot_index:2*nband_in**2+bd2tot_index)
1088 :
1089 : ! turn el_veloc to cartesian coordinates
1090 4637 : do iband=1,nband_in
1091 17256 : do l1=1,3
1092 54704 : do ii=1,3
1093 : elph_tr_ds%el_veloc(ikpt,iband,l1,isppol)=elph_tr_ds%el_veloc(ikpt,iband,l1,isppol)+&
1094 37872 : & hdr1%rprimd(l1,ii)*eig1_k(2*iband-1+(iband-1)*2*nband_in,ii)/two_pi
1095 : im_el_veloc(l1)=im_el_veloc(l1)+&
1096 12624 : & hdr1%rprimd(l1,ii)*eig1_k(2*iband+(iband-1)*2*nband_in,ii)/two_pi
1097 : end do
1098 : end do ! l1
1099 : end do
1100 : end do
1101 : end do ! end isppol
1102 :
1103 5 : call krank%free()
1104 5 : ABI_FREE(eig1_k)
1105 5 : ABI_FREE(eigen11)
1106 5 : ABI_FREE(eigen12)
1107 5 : ABI_FREE(eigen13)
1108 :
1109 5 : call hdr1%free()
1110 :
1111 5 : write(std_out,*)'out of read_el_veloc '
1112 :
1113 5 : end subroutine read_el_veloc
1114 : !!***
1115 :
1116 : !!****f* m_iogkk/inpgkk
1117 : !! NAME
1118 : !! inpgkk
1119 : !!
1120 : !! FUNCTION
1121 : !! read in gkk file and return eigenvalue matrix
1122 : !! Only works for a single gkk matrix (1 perturbation and qpoint) in the file
1123 : !! like the files produced by outgkk
1124 : !!
1125 : !! INPUTS
1126 : !!
1127 : !! filegkk= filename
1128 : !!
1129 : !! OUTPUT
1130 : !! eigen1 = response function 1st order eigenvalue matrix
1131 : !!
1132 : !! SOURCE
1133 :
1134 15 : subroutine inpgkk(eigen1,filegkk,hdr1)
1135 :
1136 : !Arguments ------------------------------------
1137 : !scalars
1138 : character(len=fnlen),intent(in) :: filegkk
1139 : type(hdr_type), intent(out) :: hdr1
1140 : !arrays
1141 : real(dp),allocatable,intent(out) :: eigen1(:)
1142 :
1143 : !Local variables-------------------------------
1144 : !scalars
1145 : integer :: bantot1
1146 : integer :: isppol, ikpt, mband, ikb
1147 : integer :: unitgkk, fform, ierr, n1wf, i1wf
1148 15 : type(hdr_type) :: hdr0
1149 15 : real(dp), allocatable :: eigen(:)
1150 : character(len=500) :: message
1151 :
1152 : ! *************************************************************************
1153 :
1154 15 : if (open_file(filegkk,message,newunit=unitgkk,form='unformatted',status='old') /= 0) then
1155 0 : ABI_ERROR(message)
1156 : end if
1157 :
1158 : !read in header of GS file and eigenvalues
1159 15 : call hdr0%fort_read(unitgkk, fform)
1160 15 : ABI_CHECK(fform /= 0, "hdr_fort_read returned fform == 0")
1161 :
1162 1287 : mband = maxval(hdr0%nband(:))
1163 45 : ABI_MALLOC(eigen,(mband))
1164 15 : call wrtout(std_out,'inpgkk : try to reread GS eigenvalues','COLL')
1165 :
1166 30 : do isppol=1,hdr0%nsppol
1167 1302 : do ikpt=1,hdr0%nkpt
1168 1272 : read (unitgkk,IOSTAT=ierr) eigen(1:hdr0%nband(ikpt))
1169 1287 : ABI_CHECK(ierr==0,'reading eigen from gkk file')
1170 : end do
1171 : end do
1172 :
1173 15 : read(unitgkk,IOSTAT=ierr) n1wf
1174 15 : ABI_CHECK(ierr==0,"reading n1wf from gkk file")
1175 :
1176 15 : ABI_FREE(eigen)
1177 15 : call hdr0%free()
1178 :
1179 15 : if (n1wf > 1) then
1180 : write(message,'(3a)')&
1181 0 : & 'several 1wf records were found in the file,',ch10, &
1182 0 : & 'which is not allowed for reading with this routine'
1183 0 : ABI_ERROR(message)
1184 : end if
1185 :
1186 : !read in header of 1WF file
1187 15 : call hdr1%fort_read(unitgkk, fform)
1188 15 : if (fform == 0) then
1189 0 : write(message,'(a,i0,a)')' 1WF header number ',i1wf,' was mis-read. fform == 0'
1190 0 : ABI_ERROR(message)
1191 : end if
1192 :
1193 15 : bantot1 = 2*hdr1%nsppol*hdr1%nkpt*mband**2
1194 45 : ABI_MALLOC(eigen1, (bantot1))
1195 :
1196 :
1197 : !retrieve 1WF <psi_k+q | H | psi_k> from gkk file and echo to output
1198 15 : ikb = 0
1199 30 : do isppol=1,hdr1%nsppol
1200 1302 : do ikpt=1,hdr1%nkpt
1201 1272 : read (unitgkk,IOSTAT=ierr) eigen1(ikb+1:ikb+2*hdr1%nband(ikpt)**2)
1202 1272 : ikb = ikb + 2*hdr1%nband(ikpt)**2
1203 1287 : if (ierr /= 0) then
1204 0 : write(message,'(a,2i0)')'reading eigen1 from gkk file, spin, kpt_idx',isppol,ikpt
1205 0 : ABI_ERROR(message)
1206 : end if
1207 : end do
1208 : end do
1209 :
1210 15 : close(unitgkk)
1211 :
1212 15 : end subroutine inpgkk
1213 : !!***
1214 :
1215 : !!****f* m_iogkk/completeperts
1216 : !!
1217 : !! NAME
1218 : !! completeperts
1219 : !!
1220 : !! FUNCTION
1221 : !! Complete perturbations wrt atoms and reduced directions
1222 : !! for a fixed qpoint. Normally there is a test in read_gkk which guarantees
1223 : !! that enough irreducible perturbations are present to generate everything.
1224 : !! h1_mat_el is first squared, making a (ipert,jpert) matrix which has the same
1225 : !! symmetry properties as the dynamical matrix.
1226 : !!
1227 : !! INPUTS
1228 : !! Cryst<crystal_t>=Info on the unit cell and symmetries.
1229 : !! nbranch=Number of phonon branches.
1230 : !! nFSband=Number of bands in H1 matrix elements.
1231 : !! nkpt=Number of k-points in matrix elements.
1232 : !! nsppol=Number of independent spin polarizations.
1233 : !! gkk_flag = flags for presence of gkk matrix elements
1234 : !! h1_mat_el = irreducible matrix elements to be completed and squared
1235 : !! qpt = qpoint
1236 : !! symq = flags for symmetry elements conserving the present qpoint
1237 : !! tnons = translation vectors associated with symops
1238 : !!
1239 : !! OUTPUT
1240 : !! h1_mat_el_sq = irreducible matrix elements squared and completed
1241 : !! gkk_flag = changed on output
1242 : !!
1243 : !! SOURCE
1244 :
1245 48 : subroutine completeperts(Cryst,nbranch,nFSband,nkpt,nsppol,gkk_flag,h1_mat_el,h1_mat_el_sq,&
1246 48 : & qpt,symq,qtimrev)
1247 :
1248 : !Arguments ------------------------------------
1249 : !scalars
1250 : integer,intent(in) :: qtimrev,nbranch,nFSband,nkpt,nsppol
1251 : type(crystal_t),intent(in) :: Cryst
1252 : !arrays
1253 : integer,intent(in) :: symq(4,2,Cryst%nsym)
1254 : integer,intent(inout) :: gkk_flag(nbranch,nbranch,nkpt,nsppol)
1255 : real(dp),intent(in) :: qpt(3)
1256 : real(dp),intent(in) :: h1_mat_el(2,nFSband**2,nbranch,nkpt,nsppol)
1257 : real(dp),intent(out) :: h1_mat_el_sq(2,nFSband**2,nbranch**2,nkpt,nsppol)
1258 :
1259 : !Local variables-------------------------------
1260 : !scalars
1261 : integer :: ikpt_phon,iatom1,iatom2,ibb,idir1,idir2,ipert1,ipert2,isppol,mpert,natom
1262 : real(dp) :: im1,im2,re1,re2,res
1263 : character(len=500) :: msg
1264 : !arrays
1265 : integer,allocatable :: tmpflg(:,:,:,:)
1266 48 : real(dp),allocatable :: tmpval(:,:,:,:,:)
1267 :
1268 : ! *************************************************************************
1269 :
1270 : !WARNING! Stupid patch in d2sym3 imposes these matrices to have size natom+2
1271 48 : natom = Cryst%natom
1272 48 : mpert = natom+2
1273 :
1274 240 : ABI_MALLOC(tmpflg,(3,mpert,3,mpert))
1275 240 : ABI_MALLOC(tmpval,(2,3,mpert,3,mpert))
1276 :
1277 7441299 : h1_mat_el_sq = zero
1278 :
1279 288 : write (std_out,*) ' completeperts: shape(h1_mat_el_sq) = ', shape(h1_mat_el_sq)
1280 :
1281 99 : do isppol=1,nsppol
1282 51 : write(std_out,*)'completeperts: isppol = ', isppol
1283 : !
1284 4203 : do ikpt_phon=1,nkpt
1285 176232 : do ibb=1,nFSband**2
1286 : !
1287 55943712 : tmpval = zero
1288 23801184 : tmpflg = 0
1289 : ! for a fixed k (q) band and sppol construct the gamma matrix for (3 natom)^2 perturbation pairs
1290 378048 : do iatom1=1,natom
1291 995808 : do idir1=1,3
1292 617760 : ipert1 = (iatom1-1)*3+idir1
1293 617760 : if (gkk_flag(ipert1,ipert1,ikpt_phon,isppol) < 0) cycle
1294 617760 : re1 = h1_mat_el(1,ibb,ipert1,ikpt_phon,isppol)
1295 617760 : im1 = h1_mat_el(2,ibb,ipert1,ikpt_phon,isppol)
1296 :
1297 1644192 : do iatom2=1,natom
1298 3899808 : do idir2=1,3
1299 2461536 : ipert2 = (iatom2-1)*3+idir2
1300 2461536 : if (gkk_flag(ipert2,ipert2,ikpt_phon,isppol) < 0) cycle
1301 2461536 : tmpflg(idir1,iatom1,idir2,iatom2) = 1
1302 2461536 : re2 = h1_mat_el(1,ibb,ipert2,ikpt_phon,isppol)
1303 2461536 : im2 = h1_mat_el(2,ibb,ipert2,ikpt_phon,isppol)
1304 : !
1305 : ! conjg(h1_mat_el_2) * h1_mat_el_1
1306 2461536 : res = re1*re2 + im1*im2
1307 2461536 : tmpval(1,idir1,iatom1,idir2,iatom2) = res
1308 2461536 : res = re1*im2 - im1*re2
1309 3282048 : tmpval(2,idir1,iatom1,idir2,iatom2) = res
1310 :
1311 : end do !idir2
1312 : end do !iatom2
1313 : end do !idir1
1314 : end do !iatom1
1315 :
1316 : ! matrix is symmetrized like a dynamical matrix. No change of band or k
1317 : ! in here. This should be checked (if we have to restrict further the symmetry operations)
1318 172128 : call d2sym3(tmpflg,tmpval,Cryst%indsym,mpert,natom,Cryst%nsym,qpt,symq,Cryst%symrec,Cryst%symrel,qtimrev,1)
1319 4277856 : if (sum(tmpflg(:,1:natom,:,1:natom)) /= 3*natom*3*natom) then
1320 : write(msg,'(3a,4i0)')&
1321 0 : & 'A perturbation is missing after completion with d2sym3',ch10,&
1322 0 : & 'tmpflg, ikpt_phon, isppol: ',tmpflg,ikpt_phon,isppol
1323 0 : ABI_ERROR(msg)
1324 : end if
1325 : !
1326 : ! Save values for calculation of |gkk|^2
1327 382152 : do iatom1=1,natom
1328 995808 : do idir1=1,3
1329 617760 : ipert1 = (iatom1-1)*3+idir1
1330 1644192 : do iatom2=1,natom
1331 3899808 : do idir2=1,3
1332 : !
1333 : ! mjv 29/10/2007 ipert2 now contains the composite index ip1*nperts+ip2
1334 2461536 : ipert2 = (iatom2-1)*3 + idir2 + (ipert1-1)*3*natom
1335 2461536 : h1_mat_el_sq(1,ibb,ipert2,ikpt_phon,isppol) = pi*tmpval(1,idir2,iatom2,idir1,iatom1)
1336 3282048 : h1_mat_el_sq(2,ibb,ipert2,ikpt_phon,isppol) = pi*tmpval(2,idir2,iatom2,idir1,iatom1)
1337 : end do
1338 : end do
1339 : end do
1340 : end do
1341 : !
1342 : end do !end ibb band dos
1343 : !
1344 : ! Set flags.
1345 18195 : do ipert1=1,3*natom
1346 70632 : do ipert2=1,3*natom
1347 66528 : if (gkk_flag(ipert2,ipert1,ikpt_phon,isppol) < 0) gkk_flag(ipert2,ipert1,ikpt_phon,isppol) = 1
1348 : end do
1349 : end do
1350 :
1351 : end do !end kpt_phon do
1352 : end do !end sppol do
1353 :
1354 48 : ABI_FREE(tmpflg)
1355 48 : ABI_FREE(tmpval)
1356 :
1357 48 : end subroutine completeperts
1358 : !!***
1359 :
1360 : !!****f* ABINIT/normsq_gkq
1361 : !!
1362 : !! NAME
1363 : !! normsq_gkq
1364 : !!
1365 : !! FUNCTION
1366 : !! This routine takes the gkq matrix elements for a given qpoint,
1367 : !! does the scalar product with the phonon displacement vector,
1368 : !! squares the gkq matrix elements multiplies by the appropriate weights
1369 : !! and puts them in a uniform (atom,icart) basis
1370 : !!
1371 : !! INPUTS
1372 : !! displ_red = phonon mode displacement vectors in reduced coordinated.
1373 : !! eigvec = eigenvectors of phonons (to turn to cartesian coord frame)
1374 : !! elph_ds = datastructure with gkk matrix elements
1375 : !! FSfullpqtofull = mapping of k + q to k
1376 : !! h1_mat_el_sq = matrix elements $<psi_{k+q,m} | H^{1} | psi_{k,n}>$ matrix-squared
1377 : !! iqptirred = index of present qpoint
1378 : !! phfrq_tmp = phonon frequencies
1379 : !! qpt_irred = array of qpoint coordinates
1380 : !!
1381 : !! OUTPUT
1382 : !! elph_ds%gkq filled
1383 : !! qdata(elph_ds%nbranch,elph_ds%nsppol,3) = array containing the phonon frequency, the linewidth and $\lambda_{q,\nu}$.
1384 : !!
1385 : !! SOURCE
1386 :
1387 48 : subroutine normsq_gkq(displ_red,eigvec,elph_ds,FSfullpqtofull,&
1388 48 : & h1_mat_el_sq,iqptirred,phfrq_tmp,qpt_irred,qdata)
1389 :
1390 : !Arguments ------------------------------------
1391 : !scalars
1392 : integer,intent(in) :: iqptirred
1393 : type(elph_type),intent(inout) :: elph_ds
1394 : !arrays
1395 : integer,intent(in) :: FSfullpqtofull(elph_ds%k_phon%nkpt,elph_ds%nqpt_full)
1396 : real(dp),intent(in) :: displ_red(2,elph_ds%nbranch,elph_ds%nbranch)
1397 : real(dp),intent(in) :: eigvec(2,elph_ds%nbranch,elph_ds%nbranch)
1398 : real(dp),intent(inout) :: &
1399 : & h1_mat_el_sq(2,elph_ds%nFSband*elph_ds%nFSband,elph_ds%nbranch*elph_ds%nbranch,elph_ds%k_phon%my_nkpt,elph_ds%nsppol)
1400 : real(dp),intent(in) :: phfrq_tmp(elph_ds%nbranch),qpt_irred(3,elph_ds%nqptirred)
1401 : real(dp),intent(out) :: qdata(elph_ds%nbranch,elph_ds%nsppol,3)
1402 :
1403 : !Local variables-------------------------------
1404 : !scalars
1405 : integer :: i1,i2,ier,ii,isppol,jbranch,comm
1406 : real(dp) :: lambda_tot
1407 : character(len=500) :: message
1408 : !arrays
1409 96 : real(dp) :: accum_mat(2,elph_ds%nbranch,elph_ds%nbranch,elph_ds%nsppol)
1410 96 : real(dp) :: accum_mat2(2,elph_ds%nbranch,elph_ds%nbranch,elph_ds%nsppol)
1411 96 : real(dp) :: gam_now2(2,elph_ds%nbranch,elph_ds%nbranch)
1412 96 : real(dp) :: lambda(elph_ds%nsppol)
1413 48 : real(dp),allocatable :: matrx(:,:),val(:),vec(:,:,:)
1414 48 : real(dp),allocatable :: zhpev1(:,:),zhpev2(:)
1415 :
1416 : ! *************************************************************************
1417 :
1418 : DBG_ENTER("COLL")
1419 :
1420 2637 : accum_mat = zero
1421 2637 : accum_mat2 = zero
1422 48 : comm = xmpi_world
1423 :
1424 48 : if (elph_ds%ep_scalprod == 1) then
1425 : !
1426 15 : if (elph_ds%ep_keepbands == 0) then
1427 0 : call wrtout(std_out,' normsq_gkq : calling nmsq_gam_sumFS',"COLL")
1428 : call nmsq_gam_sumFS (accum_mat,accum_mat2,displ_red,eigvec,elph_ds,FSfullpqtofull,&
1429 0 : & h1_mat_el_sq,iqptirred)
1430 :
1431 15 : else if (elph_ds%ep_keepbands == 1) then
1432 15 : call wrtout(std_out,' normsq_gkq : calling nmsq_gam',"COLL")
1433 : call nmsq_gam (accum_mat,accum_mat2,displ_red,eigvec,elph_ds,FSfullpqtofull,&
1434 15 : & h1_mat_el_sq,iqptirred)
1435 :
1436 : else
1437 0 : write (message,'(a,i0)')' Wrong value for elph_ds%ep_keepbands = ',elph_ds%ep_keepbands
1438 0 : ABI_BUG(message)
1439 : end if
1440 : !
1441 33 : else if (elph_ds%ep_scalprod == 0) then ! Interpolate on the pure "matrix of matrix elements" and do the scalar products later.
1442 : !
1443 33 : if (elph_ds%ep_keepbands == 0) then
1444 30 : call wrtout(std_out,' normsq_gkq : calling nmsq_pure_gkk_sumFS',"COLL")
1445 : call nmsq_pure_gkk_sumFS (accum_mat,accum_mat2,displ_red,elph_ds,FSfullpqtofull,&
1446 30 : & h1_mat_el_sq,iqptirred)
1447 :
1448 3 : else if (elph_ds%ep_keepbands == 1) then
1449 3 : call wrtout(std_out,' normsq_gkq : calling nmsq_pure_gkk',"COLL")
1450 :
1451 : call nmsq_pure_gkk (accum_mat,accum_mat2,displ_red,elph_ds,FSfullpqtofull,&
1452 3 : & h1_mat_el_sq,iqptirred)
1453 : else
1454 0 : write (message,'(a,i0)')' Wrong value for elph_ds%ep_keepbands = ',elph_ds%ep_keepbands
1455 0 : ABI_BUG(message)
1456 : end if
1457 : !
1458 : else
1459 0 : write (message,'(a,i0)')' Wrong value for elph_ds%ep_scalprod = ',elph_ds%ep_scalprod
1460 0 : ABI_BUG(message)
1461 : end if
1462 : !end if flag for doing scalar product now.
1463 :
1464 :
1465 : !MG: values without the good prefactor
1466 2637 : accum_mat = accum_mat * elph_ds%occ_factor/elph_ds%k_phon%nkpt
1467 :
1468 : !MG: accum_mat2 contains the line-widhts before the Fourier interpolation
1469 2637 : accum_mat2 = accum_mat2 * elph_ds%occ_factor/elph_ds%k_phon%nkpt
1470 :
1471 : !mpi sum over procs for accum_mat2
1472 48 : call xmpi_sum (accum_mat, comm, ier)
1473 48 : call xmpi_sum (accum_mat2, comm, ier)
1474 :
1475 : !MG20060531i
1476 : !write e-ph quantities before Fourier interpolation
1477 : !save e-ph values in the temporary array qdata that will be copied into elph_ds%qgrid_data
1478 :
1479 48 : write (message,'(4a,3es16.6,63a)')ch10, &
1480 48 : & ' Phonon linewidths before interpolation ',ch10, &
1481 2976 : & ' Q point = ',qpt_irred(:,iqptirred),ch10,('=',ii=1,60),ch10,&
1482 3024 : & ' Mode Frequency (Ha) Linewidth (Ha) Lambda '
1483 48 : call wrtout(std_out,message,'COLL')
1484 :
1485 48 : lambda_tot = zero
1486 99 : do isppol=1,elph_ds%nsppol
1487 288 : do ii=1,elph_ds%nbranch
1488 189 : lambda(isppol)=zero
1489 : ! MG: the tolerance factor is somehow arbitrary
1490 189 : if (abs(phfrq_tmp(ii)) > tol10) lambda(isppol)=accum_mat2(1,ii,ii,isppol)/&
1491 141 : & (pi*elph_ds%n0(isppol)*phfrq_tmp(ii)**2)
1492 189 : lambda_tot=lambda_tot+lambda(isppol)
1493 189 : write(message,'(i8,es20.6,2es16.6)' )ii,phfrq_tmp(ii),accum_mat2(1,ii,ii,isppol),lambda(isppol)
1494 189 : call wrtout(std_out,message,'COLL')
1495 : ! save values
1496 189 : qdata(ii,isppol,1)=phfrq_tmp(ii)
1497 189 : qdata(ii,isppol,2)=accum_mat2(1,ii,ii,isppol)
1498 240 : qdata(ii,isppol,3)=lambda(isppol)
1499 : end do !loop over branch
1500 : end do !loop over sppol
1501 :
1502 : !normalize for number of spins
1503 48 : lambda_tot = lambda_tot / elph_ds%nsppol
1504 :
1505 5808 : write(message,'(61a,44x,es16.6,62a)' )('=',ii=1,60),ch10,lambda_tot,ch10,('=',ii=1,60),ch10
1506 48 : call wrtout(std_out,message,'COLL')
1507 : !ENDMG20060531
1508 :
1509 : !immediately calculate linewidths:
1510 48 : write(std_out,*) 'summed accum_mat = '
1511 48 : write(std_out,'(3(2E18.6,1x))') accum_mat(:,:,:,1)
1512 48 : write(std_out,*) 'summed accum_mat2 = '
1513 228 : write(std_out,'(3(2E18.6,1x))') (accum_mat2(:,ii,ii,1),ii=1,elph_ds%nbranch)
1514 48 : write(std_out,*) 'displ_red = '
1515 48 : write(std_out,'(3(2E18.6,1x))') displ_red
1516 :
1517 48 : if (elph_ds%ep_scalprod == 1) then
1518 30 : do isppol=1,elph_ds%nsppol
1519 : ! Diagonalize gamma matrix at qpoint (complex matrix). Copied from dfpt_phfrq
1520 15 : ier=0
1521 15 : ii=1
1522 45 : ABI_MALLOC(matrx,(2,(elph_ds%nbranch*(elph_ds%nbranch+1))/2))
1523 96 : do i2=1,elph_ds%nbranch
1524 366 : do i1=1,i2
1525 270 : matrx(1,ii)=accum_mat2(1,i1,i2,isppol)
1526 270 : matrx(2,ii)=accum_mat2(2,i1,i2,isppol)
1527 351 : ii=ii+1
1528 : end do
1529 : end do
1530 45 : ABI_MALLOC(zhpev1,(2,2*elph_ds%nbranch-1))
1531 45 : ABI_MALLOC(zhpev2,(3*elph_ds%nbranch-2))
1532 45 : ABI_MALLOC(val,(elph_ds%nbranch))
1533 60 : ABI_MALLOC(vec,(2,elph_ds%nbranch,elph_ds%nbranch))
1534 15 : call ZHPEV ('V','U',elph_ds%nbranch,matrx,val,vec,elph_ds%nbranch,zhpev1,zhpev2,ier)
1535 :
1536 15 : write (std_out,*) ' normsq_gkq : accumulated eigenvalues isppol ',isppol, ' = '
1537 15 : write (std_out,'(3E18.6)') val
1538 15 : ABI_FREE(matrx)
1539 15 : ABI_FREE(zhpev1)
1540 15 : ABI_FREE(zhpev2)
1541 15 : ABI_FREE(vec)
1542 30 : ABI_FREE(val)
1543 : end do ! isppol
1544 :
1545 33 : else if (elph_ds%ep_scalprod == 0) then
1546 :
1547 :
1548 69 : do isppol=1,elph_ds%nsppol
1549 36 : call gam_mult_displ(elph_ds%nbranch, displ_red, accum_mat(:,:,:,isppol), gam_now2)
1550 :
1551 36 : write (std_out,*) ' normsq_gkq : accumulated eigenvalues isppol ', isppol, ' = '
1552 144 : write (std_out,'(3(E14.6,1x))') (gam_now2(1,jbranch,jbranch), jbranch=1,elph_ds%nbranch)
1553 36 : write (std_out,*) ' normsq_gkq : imag part = '
1554 177 : write (std_out,'(3(E14.6,1x))') (gam_now2(2,jbranch,jbranch), jbranch=1,elph_ds%nbranch)
1555 : end do ! isppol
1556 :
1557 : end if
1558 :
1559 : DBG_EXIT("COLL")
1560 :
1561 48 : end subroutine normsq_gkq
1562 : !!***
1563 :
1564 : !!****f* ABINIT/nmsq_gam
1565 : !!
1566 : !! NAME
1567 : !! nmsq_gam
1568 : !!
1569 : !! FUNCTION
1570 : !! Calculate gamma matrices keeping full dependence on bands
1571 : !! from original h1_mat_el_sq matrix elements (no averaging over
1572 : !! bands near the Fermi surface)
1573 : !!
1574 : !! INPUTS
1575 : !! displ_red = phonon mode displacement vectors, post-multiplied by gprim matrix
1576 : !! (ie. turned to reduced coordinates)
1577 : !! eigvec = phonon eigenvectors
1578 : !! elph_ds = datastructure with gkk matrix elements
1579 : !! FSfullpqtofull = mapping of k+q to k
1580 : !! kpt_phon = coordinates of kpoints near to FS
1581 : !! h1_mat_el_sq = matrix elements $<psi_{k+q,m} | H^{1} | psi_{k,n}>$ squared
1582 : !! iqptirred = index of present qpoint
1583 : !!
1584 : !! OUTPUT
1585 : !! accum_mat = matrix for accumulating FS average of gkk (gamma matrix -> linewidths)
1586 : !! accum_mat2 = matrix for accumulating FS average of gamma matrix with good prefactors
1587 : !!
1588 : !! SOURCE
1589 :
1590 15 : subroutine nmsq_gam (accum_mat,accum_mat2,displ_red,eigvec,elph_ds,FSfullpqtofull,&
1591 15 : & h1_mat_el_sq,iqptirred)
1592 :
1593 : !Arguments ------------------------------------
1594 : !scalars
1595 : integer,intent(in) :: iqptirred
1596 : type(elph_type),intent(inout) :: elph_ds
1597 : !arrays
1598 : integer,intent(in) :: FSfullpqtofull(elph_ds%k_phon%nkpt,elph_ds%nqpt_full)
1599 : real(dp),intent(in) :: displ_red(2,elph_ds%nbranch,elph_ds%nbranch)
1600 : real(dp),intent(in) :: eigvec(2,elph_ds%nbranch,elph_ds%nbranch)
1601 : real(dp),intent(inout) :: &
1602 : & h1_mat_el_sq(2,elph_ds%nFSband*elph_ds%nFSband,elph_ds%nbranch*elph_ds%nbranch,elph_ds%k_phon%my_nkpt,elph_ds%nsppol)
1603 : real(dp),intent(inout) :: accum_mat(2,elph_ds%nbranch,elph_ds%nbranch,elph_ds%nsppol)
1604 : real(dp),intent(inout) :: accum_mat2(2,elph_ds%nbranch,elph_ds%nbranch,elph_ds%nsppol)
1605 :
1606 : !Local variables-------------------------------
1607 : ! tmp variables for diagonalization
1608 : !scalars
1609 : integer :: ikpt_phon,ikpt_phonq,ib1,ib2,ibeff,ibranch,isppol,ipert1
1610 : integer :: jbranch
1611 : integer :: iqpt_fullbz
1612 : integer :: ik_this_proc
1613 : real(dp) :: sd1,sd2
1614 : character(len=500) :: message
1615 : !arrays
1616 30 : real(dp) :: gkq_1band(2,elph_ds%nbranch,elph_ds%nbranch)
1617 30 : real(dp) :: tmp_mat2(2,elph_ds%nbranch,elph_ds%nbranch)
1618 30 : real(dp) :: zgemm_tmp_mat(2,elph_ds%nbranch,elph_ds%nbranch)
1619 :
1620 : ! *************************************************************************
1621 :
1622 15 : if (elph_ds%ep_keepbands == 0) then
1623 0 : write (message,'(a,i0)')' elph_ds%ep_keepbands should be 1 while is ',elph_ds%ep_keepbands
1624 0 : ABI_ERROR(message)
1625 : end if
1626 :
1627 : !MG20060603 NOTE:
1628 : !accum_mat and accum_mat2 are real, the imaginary part is used for debugging purpose
1629 : !accum_mat2 is used to store the phonon-linewidhts before interpolation
1630 :
1631 15 : iqpt_fullbz = elph_ds%qirredtofull(iqptirred)
1632 15 : write(std_out,*) 'nmsq_gam : iqptirred = ', iqptirred
1633 :
1634 30 : do isppol=1,elph_ds%nsppol
1635 798 : do ik_this_proc =1, elph_ds%k_phon%my_nkpt
1636 768 : ikpt_phon = elph_ds%k_phon%my_ikpt(ik_this_proc)
1637 :
1638 768 : ikpt_phonq = FSfullpqtofull(ikpt_phon,iqpt_fullbz)
1639 :
1640 6287 : do ib1=1,elph_ds%nFSband
1641 5504 : sd1 = elph_ds%k_phon%wtk(ib1,ikpt_phon,isppol) !weights for distance from the fermi surface
1642 :
1643 46976 : do ib2=1,elph_ds%nFSband
1644 40704 : sd2 = elph_ds%k_phon%wtk(ib2,ikpt_phonq,isppol) !weights for distance from the fermi surface
1645 40704 : ibeff = ib2+elph_ds%nFSband*(ib1-1)
1646 :
1647 4100352 : gkq_1band(:,:,:) = zero
1648 :
1649 162816 : zgemm_tmp_mat= reshape (h1_mat_el_sq(:,ibeff,:,ik_this_proc,isppol),(/2,elph_ds%nbranch,elph_ds%nbranch/))
1650 :
1651 40704 : call gam_mult_displ(elph_ds%nbranch, displ_red, zgemm_tmp_mat, tmp_mat2)
1652 :
1653 : ! sum over bands
1654 264192 : do ipert1=1,elph_ds%nbranch
1655 264192 : gkq_1band(1,ipert1,ipert1) = gkq_1band(1,ipert1,ipert1) + tmp_mat2(1,ipert1,ipert1)
1656 : end do
1657 :
1658 : ! summing over k points and bands, still diagonal in jbranch
1659 4100352 : accum_mat(:,:,:,isppol) = accum_mat(:,:,:,isppol) + gkq_1band(:,:,:)*sd1*sd2
1660 :
1661 : ! MG20060603 : summing over bands and kpoints with weights to calculate the phonon linewidth
1662 264192 : do jbranch=1,elph_ds%nbranch
1663 711168 : accum_mat2(:,jbranch,jbranch,isppol) = accum_mat2(:,jbranch,jbranch,isppol) + gkq_1band(:,jbranch,jbranch)*sd1*sd2
1664 : end do
1665 : ! END MG
1666 :
1667 :
1668 : ! now turn to cartesian coordinates
1669 :
1670 : ! Final Gamma matrix (hermitian) = E * D_g * E^{+}
1671 : ! Where E^{+} is the hermitian conjugate of the eigenvector matrix E
1672 : ! And D_g is the diagonal matrix of values of gamma for this qpoint
1673 :
1674 : ! Here gkq_1band is indexed with real phonon modes (not atom+idir)
1675 : ! turn gkq_1band to atom+cartesian coordinates (instead of normal coordinates for qpoint)
1676 4100352 : tmp_mat2(:,:,:) = zero
1677 264192 : do ibranch =1,elph_ds%nbranch
1678 1542912 : do jbranch =1,elph_ds%nbranch
1679 : tmp_mat2(1,ibranch,jbranch) = tmp_mat2(1,ibranch,jbranch) + &
1680 1278720 : & eigvec(1,ibranch,jbranch) * gkq_1band(1,jbranch,jbranch)
1681 : tmp_mat2(2,ibranch,jbranch) = tmp_mat2(2,ibranch,jbranch) + &
1682 1502208 : & eigvec(2,ibranch,jbranch) * gkq_1band(1,jbranch,jbranch)
1683 : end do
1684 : end do
1685 4100352 : gkq_1band(:,:,:) = zero
1686 :
1687 : ! here eigvec is transposed and complexconjugated.
1688 4100352 : zgemm_tmp_mat=zero
1689 : call zgemm('n','c',elph_ds%nbranch,elph_ds%nbranch,elph_ds%nbranch,cone,&
1690 40704 : & tmp_mat2,elph_ds%nbranch,eigvec,elph_ds%nbranch,czero,zgemm_tmp_mat,elph_ds%nbranch)
1691 :
1692 4100352 : gkq_1band = zgemm_tmp_mat
1693 :
1694 : ! gamma matrix contribution in cartesian coordinates (ie interpolatable form)
1695 3963776 : h1_mat_el_sq(:,ibeff,:,ik_this_proc,isppol) = reshape(gkq_1band,(/2,elph_ds%nbranch*elph_ds%nbranch/))
1696 :
1697 : end do
1698 : end do
1699 : ! END loop over bands ib1 ib2
1700 :
1701 : end do
1702 : ! END loop over kpt_phon
1703 : end do
1704 : !END loop over nsppol
1705 :
1706 :
1707 15 : end subroutine nmsq_gam
1708 : !!***
1709 :
1710 : !!****f* ABINIT/nmsq_gam_sumfs
1711 : !!
1712 : !! NAME
1713 : !! nmsq_gam_sumfs
1714 : !!
1715 : !! FUNCTION
1716 : !! Calculate gamma matrices from original h1_mat_el_sq matrix
1717 : !! elements averaging over bands near the Fermi surface
1718 : !!
1719 : !! INPUTS
1720 : !! displ_red = phonon mode displacement vectors, post-multiplied by gprim matrix
1721 : !! (ie. turned to reduced coordinates)
1722 : !! eigvec = eigenvectors of phonons (to turn to cartesian coord frame)
1723 : !! elph_ds = datastructure with gkk matrix elements
1724 : !! FSfullpqtofull = mapping of k+q to k
1725 : !! kpt_phon = coordinates of kpoints near to FS
1726 : !! h1_mat_el_sq = matrix elements $<psi_{k+q,m} | H^{1} | psi_{k,n}>$ squared
1727 : !! iqptirred = index of present qpoint
1728 : !!
1729 : !! OUTPUT
1730 : !! accum_mat = matrix for accumulating FS average of gkk (gamma matrix -> linewidths)
1731 : !! accum_mat2 = matrix for accumulating FS average of gamma matrix with good prefactors
1732 : !!
1733 : !! SOURCE
1734 :
1735 0 : subroutine nmsq_gam_sumFS(accum_mat,accum_mat2,displ_red,eigvec,elph_ds,FSfullpqtofull,&
1736 0 : & h1_mat_el_sq,iqptirred)
1737 :
1738 : !Arguments ------------------------------------
1739 : !scalars
1740 : integer,intent(in) :: iqptirred
1741 : type(elph_type),intent(inout) :: elph_ds
1742 : !arrays
1743 : integer,intent(in) :: FSfullpqtofull(elph_ds%k_phon%nkpt,elph_ds%nqpt_full)
1744 : real(dp),intent(in) :: displ_red(2,elph_ds%nbranch,elph_ds%nbranch)
1745 : real(dp),intent(in) :: eigvec(2,elph_ds%nbranch,elph_ds%nbranch)
1746 : real(dp),intent(inout) :: &
1747 : & h1_mat_el_sq(2,elph_ds%nFSband*elph_ds%nFSband,elph_ds%nbranch*elph_ds%nbranch,elph_ds%k_phon%my_nkpt,elph_ds%nsppol)
1748 : real(dp),intent(inout) :: accum_mat(2,elph_ds%nbranch,elph_ds%nbranch,elph_ds%nsppol)
1749 : real(dp),intent(inout) :: accum_mat2(2,elph_ds%nbranch,elph_ds%nbranch,elph_ds%nsppol)
1750 :
1751 : !Local variables-------------------------------
1752 : !scalars
1753 : integer :: ikpt_phon,ikpt_phonq,ib1,ib2,ibeff,ibranch,ipert1,isppol,jbranch,iqpt_fullbz
1754 : integer :: ik_this_proc
1755 : real(dp) :: sd1,sd2
1756 : character(len=500) :: message
1757 : !arrays
1758 0 : real(dp) :: gkq_sum_bands(2,elph_ds%nbranch,elph_ds%nbranch)
1759 0 : real(dp) :: tmp_gkq_sum_bands(2,elph_ds%nbranch,elph_ds%nbranch)
1760 0 : real(dp) :: tmp_mat2(2,elph_ds%nbranch,elph_ds%nbranch)
1761 0 : real(dp),allocatable :: zgemm_tmp_mat(:,:,:)
1762 :
1763 : ! *************************************************************************
1764 :
1765 0 : if (elph_ds%ep_keepbands /= 0) then
1766 0 : write (message,'(a,i0)')' elph_ds%ep_keepbands should be 0 in order to average over bands!',elph_ds%ep_keepbands
1767 0 : ABI_ERROR(message)
1768 : end if
1769 :
1770 0 : iqpt_fullbz = elph_ds%qirredtofull(iqptirred)
1771 :
1772 :
1773 : !MG20060603 NOTE:
1774 : !accum_mat and accum_mat2 are real, the imaginary part is used for debugging purpose
1775 : !accum_mat2 is used to store the phonon-linewidhts before interpolation
1776 :
1777 0 : ABI_MALLOC(zgemm_tmp_mat ,(2,elph_ds%nbranch,elph_ds%nbranch))
1778 :
1779 0 : do isppol=1,elph_ds%nsppol
1780 0 : do ik_this_proc =1, elph_ds%k_phon%my_nkpt
1781 0 : ikpt_phon = elph_ds%k_phon%my_ikpt(ik_this_proc)
1782 :
1783 0 : ikpt_phonq = FSfullpqtofull(ikpt_phon,iqpt_fullbz)
1784 :
1785 0 : gkq_sum_bands = zero
1786 0 : tmp_gkq_sum_bands = zero
1787 :
1788 :
1789 0 : do ib1=1,elph_ds%nFSband
1790 : ! weights for distance from the fermi surface
1791 0 : sd1 = elph_ds%k_phon%wtk(ib1,ikpt_phon,isppol)
1792 :
1793 0 : do ib2=1,elph_ds%nFSband
1794 : ! weights for distance from the fermi surface
1795 0 : sd2 = elph_ds%k_phon%wtk(ib2,ikpt_phonq,isppol)
1796 0 : ibeff=ib2+(ib1-1)*elph_ds%nFSband
1797 :
1798 0 : zgemm_tmp_mat = reshape(h1_mat_el_sq(:,ibeff,:,ik_this_proc,isppol),(/2,elph_ds%nbranch,elph_ds%nbranch/))
1799 :
1800 0 : call gam_mult_displ(elph_ds%nbranch, displ_red, zgemm_tmp_mat, tmp_mat2)
1801 :
1802 : ! sum over bands in gkq_sum_bands
1803 0 : do ipert1=1,elph_ds%nbranch
1804 0 : gkq_sum_bands(1,ipert1,ipert1) = gkq_sum_bands(1,ipert1,ipert1) + sd1*sd2*tmp_mat2(1,ipert1,ipert1)
1805 : end do
1806 :
1807 :
1808 :
1809 : end do
1810 : end do
1811 : ! END loop over bands
1812 :
1813 : ! summing over k points, still diagonal in jbranch
1814 0 : accum_mat(:,:,:,isppol) = accum_mat(:,:,:,isppol) + gkq_sum_bands(:,:,:)
1815 0 : accum_mat2(:,:,:,isppol) = accum_mat2(:,:,:,isppol) + gkq_sum_bands(:,:,:)
1816 :
1817 : ! summed over bands, now turn to cartesian coordinates
1818 :
1819 : ! Final Gamma matrix (hermitian) = E * D_g * E^{+}
1820 : ! Where E^{+} is the hermitian conjugate of the eigenvector matrix E
1821 : ! And D_g is the diagonal matrix of values of gamma for this qpoint
1822 :
1823 : ! Here gkq_sum_bands is indexed with real phonon modes (not atom+idir)
1824 : ! turn gkq_sum_bands to atom+cartesian coordinates (instead of normal coordinates for qpoint)
1825 : ! This is not a full matrix multiplication, just vector one, by
1826 : ! gkq_sum_bands(1,jbranch,jbranch)
1827 0 : tmp_mat2(:,:,:) = zero
1828 0 : do ibranch =1,elph_ds%nbranch
1829 0 : do jbranch =1,elph_ds%nbranch
1830 : tmp_mat2(1,ibranch,jbranch) = tmp_mat2(1,ibranch,jbranch) + &
1831 : & eigvec(1,ibranch,jbranch) * &
1832 0 : & gkq_sum_bands(1,jbranch,jbranch)
1833 : tmp_mat2(2,ibranch,jbranch) = tmp_mat2(2,ibranch,jbranch) + &
1834 : & eigvec(2,ibranch,jbranch) * &
1835 0 : & gkq_sum_bands(1,jbranch,jbranch)
1836 : end do
1837 : end do
1838 :
1839 : ! here eigvec is transposed and complex conjugated.
1840 0 : zgemm_tmp_mat=zero
1841 : call zgemm('n','c',elph_ds%nbranch,elph_ds%nbranch,elph_ds%nbranch,cone,&
1842 0 : & tmp_mat2,elph_ds%nbranch,eigvec,elph_ds%nbranch,czero,zgemm_tmp_mat,elph_ds%nbranch)
1843 :
1844 0 : gkq_sum_bands = zgemm_tmp_mat
1845 :
1846 : ! ! gamma matrix contribution in cartesian coordinates (ie interpolatable form)
1847 : ! gamma matrix contribution in reduced coordinates (ie interpolatable form)
1848 0 : h1_mat_el_sq(:,1,:,ik_this_proc,isppol) = reshape(gkq_sum_bands(:,:,:),(/2,elph_ds%nbranch*elph_ds%nbranch/))
1849 :
1850 : ! accum_mat(:,:,:,isppol) = accum_mat(:,:,:,isppol) + gkq_sum_bands(:,:,:)
1851 : end do
1852 : ! END loop over kpt_phon
1853 : end do
1854 : !END loop over sppol
1855 :
1856 0 : ABI_FREE(zgemm_tmp_mat)
1857 :
1858 0 : end subroutine nmsq_gam_sumFS
1859 : !!***
1860 :
1861 :
1862 : !!****f* ABINIT/nmsq_pure_gkk
1863 : !!
1864 : !! NAME
1865 : !! nmsq_pure_gkk
1866 : !!
1867 : !! FUNCTION
1868 : !! Calculate gamma matrices for pure gkk case, ie when the
1869 : !! scalar product with the displacement vector is done later
1870 : !! Sum over bands is carried out later.
1871 : !!
1872 : !! INPUTS
1873 : !! displ_red = phonon displacement in reduced coordinates (used to calculate the ph linewidth)
1874 : !! elph_ds = datastructure with gkk matrix elements
1875 : !! FSfullpqtofull = mapping of k+q to k
1876 : !! kpt_phon = coordinates of kpoints near to FS
1877 : !! h1_mat_el_sq = matrix elements $<psi_{k+q,m} | H^{1} | psi_{k,n}>$ squared
1878 : !! iqptirred = index of present qpoint
1879 : !!
1880 : !! OUTPUT
1881 : !! elph_ds%gkq filled
1882 : !! accum_mat = matrix for accumulating FS average of gkk (gamma matrix -> linewidths)
1883 : !! accum_mat2 = complex array whose real part contains the phonon linewidth
1884 : !!
1885 : !! SOURCE
1886 :
1887 3 : subroutine nmsq_pure_gkk(accum_mat,accum_mat2,displ_red,elph_ds,FSfullpqtofull,&
1888 3 : & h1_mat_el_sq,iqptirred)
1889 :
1890 : !Arguments ------------------------------------
1891 : !scalars
1892 : integer,intent(in) :: iqptirred
1893 : type(elph_type),intent(inout) :: elph_ds
1894 : !arrays
1895 : integer,intent(in) :: FSfullpqtofull(elph_ds%k_phon%nkpt,elph_ds%nqpt_full)
1896 : real(dp),intent(in) :: displ_red(2,elph_ds%nbranch,elph_ds%nbranch)
1897 : real(dp),intent(inout) :: &
1898 : & h1_mat_el_sq(2,elph_ds%nFSband*elph_ds%nFSband,elph_ds%nbranch*elph_ds%nbranch,elph_ds%k_phon%my_nkpt,elph_ds%nsppol)
1899 : real(dp),intent(inout) :: accum_mat(2,elph_ds%nbranch,elph_ds%nbranch,elph_ds%nsppol)
1900 : real(dp),intent(inout) :: accum_mat2(2,elph_ds%nbranch,elph_ds%nbranch,elph_ds%nsppol)
1901 :
1902 : !Local variables-------------------------------
1903 : !scalars
1904 : integer :: ikpt_phon,ikpt_phonq,ib1,ib2,ibeff,ipert1,isppol
1905 : integer :: iqpt_fullbz
1906 : integer :: ik_this_proc
1907 : real(dp) :: sd1,sd2
1908 : character(len=500) :: message
1909 : !arrays
1910 6 : real(dp) :: gkq_sum_bands(2,elph_ds%nbranch,elph_ds%nbranch)
1911 6 : real(dp) :: tmp_mat2(2,elph_ds%nbranch,elph_ds%nbranch)
1912 6 : real(dp) :: zgemm_tmp_mat(2,elph_ds%nbranch,elph_ds%nbranch)
1913 :
1914 : ! *************************************************************************
1915 :
1916 3 : if (elph_ds%ep_keepbands /= 1) then
1917 0 : message = ' elph_ds%ep_keepbands should be 1 to keep bands!'
1918 0 : ABI_ERROR(message)
1919 : end if
1920 :
1921 3 : iqpt_fullbz = elph_ds%qirredtofull(iqptirred)
1922 :
1923 : !h1_mat_el_sq is already fine here - nothing to do
1924 :
1925 :
1926 : !MG20060603 NOTE:
1927 : !accum_mat and accum_mat2 are real, the imaginary part is used for debugging purpose
1928 : !accum_mat2 is used to store the phonon-linewidhts before interpolation
1929 :
1930 : !MJV 20070525 NOTE:
1931 : !in some of the nmsq routines, in particular this one, the work done to
1932 : !calculate accum_mat,accum_mat2 is completely superfluous and will be re-done
1933 : !on the interpolated values.
1934 : !MG uses them for the QPT output, however, so keep it for consistency for the
1935 : !moment.
1936 :
1937 6 : do isppol=1,elph_ds%nsppol
1938 651 : do ik_this_proc =1, elph_ds%k_phon%my_nkpt
1939 648 : ikpt_phon = elph_ds%k_phon%my_ikpt(ik_this_proc)
1940 :
1941 648 : ikpt_phonq = FSfullpqtofull(ikpt_phon,iqpt_fullbz)
1942 :
1943 20088 : gkq_sum_bands(:,:,:) = zero
1944 :
1945 : ! gkq_sum_bands = \sum_{ib1,ib2} \langle k+q \mid H^{(1)}_{q,\tau_i,\alpha_i} \mid k \rangle
1946 : ! \cdot \langle k \mid H^{(1)}_{q,\tau_j,\alpha_j} \mid k+q \rangle
1947 : ! where ibranch -> \tau_i,\alpha_i and jbranch -> \tau_j,\alpha_j
1948 :
1949 4536 : do ib1=1,elph_ds%nFSband
1950 :
1951 3888 : sd1 = elph_ds%k_phon%wtk(ib1,ikpt_phon,isppol) ! weights for distance from the fermi surface
1952 :
1953 27864 : do ib2=1,elph_ds%nFSband
1954 :
1955 23328 : sd2 = elph_ds%k_phon%wtk(ib2,ikpt_phonq,isppol) ! weights for distance from the fermi surface
1956 23328 : ibeff = ib2+(ib1-1)*elph_ds%nFSband
1957 :
1958 : gkq_sum_bands = gkq_sum_bands + &
1959 797040 : & sd1*sd2*reshape(h1_mat_el_sq(:,ibeff,:,ik_this_proc,isppol),(/2,elph_ds%nbranch,elph_ds%nbranch/))
1960 :
1961 : end do !ib2
1962 : end do !ib1
1963 : ! END loops over bands
1964 :
1965 :
1966 20091 : accum_mat(:,:,:,isppol) = accum_mat(:,:,:,isppol) + gkq_sum_bands(:,:,:)
1967 : end do
1968 : ! END loop over kpt_phon
1969 :
1970 : ! MG20060603
1971 : ! do scalar product with the displ_red to calculate the ph lwdth before interpolation (stored in accum_mat2)
1972 :
1973 93 : zgemm_tmp_mat = accum_mat(:,:,:,isppol)
1974 :
1975 3 : call gam_mult_displ(elph_ds%nbranch, displ_red, zgemm_tmp_mat, tmp_mat2)
1976 :
1977 15 : do ipert1=1,elph_ds%nbranch
1978 12 : accum_mat2(1,ipert1,ipert1,isppol) = accum_mat2(1,ipert1,ipert1,isppol) + tmp_mat2(1,ipert1,ipert1)
1979 : end do
1980 :
1981 : ! ENDMG
1982 :
1983 : end do ! isppol
1984 :
1985 3 : end subroutine nmsq_pure_gkk
1986 : !!***
1987 :
1988 : !!****f* ABINIT/nmsq_pure_gkk_sumfs
1989 : !!
1990 : !! NAME
1991 : !! nmsq_pure_gkk_sumfs
1992 : !!
1993 : !! FUNCTION
1994 : !! Calculate gamma matrices for pure gkk case, i.e, when the
1995 : !! scalar product with the displacement vector is done later
1996 : !! Sum over bands is carried out now.
1997 : !!
1998 : !! INPUTS
1999 : !! displ_red = phonon displacement in reduced coordinates (used to calculate the ph linewidth)
2000 : !! elph_ds = datastructure with gkk matrix elements
2001 : !! FSfullpqtofull = mapping of k+q to k
2002 : !! kpt_phon = coordinates of kpoints near to FS
2003 : !! h1_mat_el_sq = matrix elements $<psi_{k+q,m} | H^{1} | psi_{k,n}>$ squared
2004 : !! iqptirred = index of present qpoint
2005 : !!
2006 : !! OUTPUT
2007 : !! accum_mat = matrix for accumulating FS average of gkk (gamma matrix -> linewidths)
2008 : !! accum_mat2 = complex array whose real part contains the phonon linewidth
2009 : !!
2010 : !! SOURCE
2011 :
2012 30 : subroutine nmsq_pure_gkk_sumfs(accum_mat,accum_mat2,displ_red,elph_ds,FSfullpqtofull,h1_mat_el_sq,iqptirred)
2013 :
2014 : !Arguments ------------------------------------
2015 : !scalars
2016 : integer,intent(in) :: iqptirred
2017 : type(elph_type),intent(in) :: elph_ds
2018 : !arrays
2019 : integer,intent(in) :: FSfullpqtofull(elph_ds%k_phon%nkpt,elph_ds%nqpt_full)
2020 : real(dp),intent(in) :: displ_red(2,elph_ds%nbranch,elph_ds%nbranch)
2021 : real(dp),intent(inout) :: &
2022 : & h1_mat_el_sq(2,elph_ds%nFSband*elph_ds%nFSband,elph_ds%nbranch*elph_ds%nbranch,elph_ds%k_phon%my_nkpt,elph_ds%nsppol)
2023 : real(dp),intent(inout) :: accum_mat(2,elph_ds%nbranch,elph_ds%nbranch,elph_ds%nsppol)
2024 : real(dp),intent(inout) :: accum_mat2(2,elph_ds%nbranch,elph_ds%nbranch,elph_ds%nsppol)
2025 :
2026 : !Local variables-------------------------------
2027 : !scalars
2028 : integer :: ikpt_phon,ikpt_phonq,ib1,ib2,ibeff,ipert1,isppol,iqpt_fullbz
2029 : integer :: nbranch,nsppol,nFSband,nkpt_phon
2030 : integer :: ik_this_proc
2031 : real(dp) :: sd1,sd2
2032 : !character(len=500) :: message
2033 : !arrays
2034 60 : real(dp) :: gkq_sum_bands(2,elph_ds%nbranch,elph_ds%nbranch)
2035 60 : real(dp) :: tmp_mat2(2,elph_ds%nbranch,elph_ds%nbranch)
2036 60 : real(dp) :: zgemm_tmp_mat(2,elph_ds%nbranch,elph_ds%nbranch)
2037 :
2038 : ! *************************************************************************
2039 :
2040 30 : if (elph_ds%ep_keepbands /= 0) then
2041 0 : ABI_BUG('ep_keepbands should be 0 to average over bands!')
2042 : end if
2043 :
2044 30 : nbranch = elph_ds%nbranch
2045 30 : nsppol = elph_ds%nsppol
2046 30 : nFSband = elph_ds%nFSband
2047 30 : nkpt_phon = elph_ds%k_phon%nkpt
2048 :
2049 30 : iqpt_fullbz = elph_ds%qirredtofull(iqptirred)
2050 :
2051 : !MG20060603 NOTE:
2052 : !accum_mat and accum_mat2 are real, the imaginary part is used for debugging purpose
2053 : !accum_mat2 is used to store the phonon-linewidhts before interpolation
2054 :
2055 63 : do isppol=1,nsppol
2056 2751 : do ik_this_proc =1, elph_ds%k_phon%my_nkpt
2057 2688 : ikpt_phon = elph_ds%k_phon%my_ikpt(ik_this_proc)
2058 :
2059 : !
2060 : ! The index of k+q in the BZ.
2061 2688 : ikpt_phonq = FSfullpqtofull(ikpt_phon,iqpt_fullbz)
2062 : !
2063 : ! gkq_sum_bands =
2064 : ! \sum_{ib1,ib2} <k+q| H^{(1)}_{q,\tau_i,\alpha_i} |k> \cdot <k| H^{(1)}_{q,\tau_j,\alpha_j}|k+q>
2065 : !
2066 : ! where ibranch = (\tau_i,\alpha_i) and jbranch = (\tau_j,\alpha_j).
2067 83328 : gkq_sum_bands(:,:,:) = zero
2068 :
2069 19488 : do ib1=1,nFSband
2070 16800 : sd1 = elph_ds%k_phon%wtk(ib1,ikpt_phon,isppol) ! weights for distance from the fermi surface
2071 :
2072 127584 : do ib2=1,nFSband
2073 108096 : sd2 = elph_ds%k_phon%wtk(ib2,ikpt_phonq,isppol) ! weights for distance from the fermi surface
2074 108096 : ibeff=ib2+(ib1-1)*nFSband
2075 :
2076 : gkq_sum_bands = gkq_sum_bands + &
2077 3692064 : & sd1*sd2* reshape(h1_mat_el_sq(:,ibeff,:,ik_this_proc,isppol),(/2,nbranch,nbranch/))
2078 : end do !ib2
2079 : end do !ib1
2080 : !
2081 : ! gamma matrix contribution in reduced coordinates (ie interpolatable form)
2082 : ! The sum over Fermi surface bands is done here, and fed into (ib1,ib2)=(1,1)
2083 80640 : h1_mat_el_sq(:,1,:,ik_this_proc,isppol) = reshape(gkq_sum_bands,(/2,nbranch**2/))
2084 :
2085 83361 : accum_mat(:,:,:,isppol) = accum_mat(:,:,:,isppol) + gkq_sum_bands(:,:,:)
2086 : end do ! kpt_phon
2087 : end do ! isppol
2088 : !
2089 : !MG20060603
2090 : !do scalar product wit displ_red to calculate the ph lwdth before interpolation (stored in accum_mat2)
2091 63 : do isppol=1,nsppol
2092 1023 : zgemm_tmp_mat = accum_mat(:,:,:,isppol)
2093 : !
2094 33 : call gam_mult_displ(nbranch, displ_red, zgemm_tmp_mat, tmp_mat2)
2095 :
2096 162 : do ipert1=1,nbranch
2097 132 : accum_mat2(1,ipert1,ipert1,isppol) = accum_mat2(1,ipert1,ipert1,isppol) + tmp_mat2(1,ipert1,ipert1)
2098 : end do
2099 : !
2100 : end do
2101 :
2102 30 : end subroutine nmsq_pure_gkk_sumfs
2103 : !!***
2104 :
2105 : end module m_iogkk
2106 : !!***
|
