Branch data 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 : : !!***
|
