Line data Source code
1 : !!****m* ABINIT/defs_elphon
2 : !!
3 : !! NAME
4 : !! defs_elphon
5 : !!
6 : !! FUNCTION
7 : !! This module contains the datastructures for elphon
8 : !! the different (huge) matrices will either be allocated and
9 : !! used, or be written to disk. All combinations should be feasible.
10 : !!
11 : !! COPYRIGHT
12 : !! Copyright (C) 2004-2026 ABINIT group (MVer, MG)
13 : !! This file is distributed under the terms of the
14 : !! GNU General Public Licence, see ~abinit/COPYING
15 : !! or http://www.gnu.org/copyleft/gpl.txt .
16 : !! For the initials of contributors, see ~abinit/doc/developers/contributors.txt .
17 : !!
18 : !! NOTES
19 : !! Contains the following datastructures:
20 : !! 1) elph_type contains data and dimensions for the kpoints near the
21 : !! fermi surface and the $g_{k k+q}$ matrix elements
22 : !!
23 : !! SOURCE
24 :
25 : #if defined HAVE_CONFIG_H
26 : #include "config.h"
27 : #endif
28 :
29 : #include "abi_common.h"
30 :
31 : module defs_elphon
32 :
33 : use defs_basis
34 : use m_abicore
35 : use m_errors
36 : use m_xmpi
37 : use m_krank
38 : use m_crystal
39 :
40 : implicit none
41 :
42 : private
43 : !!***
44 :
45 : public :: gam_mult_displ
46 : public :: complete_gamma
47 : public :: complete_gamma_tr
48 : public :: mkqptequiv
49 :
50 : !----------------------------------------------------------------------
51 : !!****t* defs_elphon/elph_kgrid_type
52 : !! NAME
53 : !! elph_kgrid_type
54 : !!
55 : !! FUNCTION
56 : !! elph_kgrid_type contains k-point grid data and dimensions
57 : !! this is a sub object of elph_type
58 : !!
59 : !! SOURCE
60 :
61 : type,public :: elph_kgrid_type
62 :
63 : integer :: nband ! number of bands for weights
64 : integer :: nsppol ! number of spin pol for weights
65 : integer :: nsym ! number of symmetry operations
66 : integer :: nkpt ! number of k-points in full grid
67 : integer :: nkptirr ! number of k-points in irreducible grid
68 : integer :: new_nkptirr ! number of k-points in irreducible grid
69 : integer :: my_nkpt ! number of k-points on present processor
70 :
71 : type(krank_t) :: krank ! ranking of all kpoints on phonon calculation grid, and inverse rank
72 :
73 : integer, allocatable :: irr2full(:) ! correspondence of irred kpoints to a full one
74 : integer, allocatable :: full2irr(:,:) ! correspondence of full k to one irred kpoints through sym and timrev
75 : integer, allocatable :: full2full(:,:,:) ! symmetry mapping of kpoints
76 : integer, allocatable :: my_kpt(:) ! flag for k-points belonging to present proc (= me index of proc for each k-point)
77 : integer, allocatable :: my_ikpt(:) ! flag for k-points belonging to present proc (= me index of proc for each k-point)
78 : integer, allocatable :: irredtoGS(:) ! (nkptirr)
79 : integer, allocatable :: new_irredtoGS(:) ! (new_nkptirr)
80 :
81 : real(dp),allocatable :: kpt(:,:) ! coordinates of the full kpoints from phonon calculation
82 : real(dp),allocatable :: kptirr(:,:) ! irreducible k-points, for preliminary set up
83 : real(dp),allocatable :: new_kptirr(:,:) ! irreducible k-points, for preliminary set up
84 : real(dp),allocatable :: wtk(:,:,:) ! integration weights (see also gkk_intweight)
85 : real(dp),allocatable :: wtq(:,:,:) ! integration weights (see also gkk_intweight)
86 : real(dp),allocatable :: wtkirr(:) ! weights for irreducible kpoints, to sum over _whole_ BZ (not just Fermi Surface)
87 : real(dp),allocatable :: new_wtkirr(:) ! weights for irreducible kpoints, to sum over _whole_ BZ (not just Fermi Surface)
88 : real(dp),allocatable :: velocwtk(:,:,:,:) ! (nFSband,nkpt_fine,3,nsppol), v(k)*wtk
89 : real(dp),allocatable :: vvelocwtk(:,:,:,:,:) ! (nFSband,nkpt_fine,3,3,nsppol), v(k)*v(k)*wtk
90 :
91 : end type elph_kgrid_type
92 :
93 : public :: elph_k_copy
94 : public :: elph_k_procs
95 : public :: elph_k_destroy
96 : !!***
97 :
98 : !----------------------------------------------------------------------
99 :
100 : !!****t* defs_elphon/elph_type
101 : !! NAME
102 : !! elph_type
103 : !!
104 : !! FUNCTION
105 : !! elph_type contains data and dimensions for the kpoints near the
106 : !! fermi surface and the $g_{k k+q}$ matrix elements
107 : !!
108 : !! SOURCE
109 :
110 : type,public :: elph_type
111 :
112 : type(elph_kgrid_type) :: k_phon ! object for k-grid of phonon calculation
113 : type(elph_kgrid_type) :: k_fine ! object for fine k-grid for FS integration
114 :
115 : integer :: natom,nbranch,nFSband,nband
116 : integer :: minFSband,maxFSband !Index of lower and upper bands used for the FS integration
117 :
118 : integer :: ngkkband !Number of bands kept in final gkk matrix elements:
119 : !either 1 if sum is performed immediately
120 : !or = nband if all elements are kept based on flag ep_keepbands
121 :
122 : integer :: nenergy ! number of points in energy
123 : ! space for the electron band energies
124 : ! energies from
125 : ! ef-nenergy*delta_e to
126 : ! ef+nenergy*delta_e
127 :
128 : integer :: n_pair ! number of pairs considered
129 :
130 : integer :: nqpt_full !number of q in full BZ
131 : integer :: nqptirred !number of irred q-points
132 :
133 :
134 : integer :: unita2f,unit_gkk2,unit_gkk_rpt
135 : integer :: unitgkq !units for file output
136 :
137 : integer :: gkqwrite
138 : integer :: gkk2write
139 : integer :: gkk_rptwrite
140 :
141 : integer :: ep_scalprod !flag to perform the scalar product
142 : integer :: symgkq !flag to symmetrize gkq matrix elements
143 : integer :: ep_keepbands !flag to sum over bands or not
144 : integer :: ep_lova ! 1 for lova, and 0 for general
145 : integer :: tuniformgrid !flag to expect uniform grid of q or not
146 : integer :: prtbltztrp !flag to output BoltzTraP input files
147 :
148 : integer :: na2f !dimensions and increments for a2F function
149 : integer :: nsppol ! number of spin polarization channels
150 : integer :: nspinor ! number of spinorial components
151 : integer :: telphint ! flag for integration over the FS with 0=tetrahedra 1=gaussians
152 : integer :: ep_nspline ! scale factor for spline interpolation in RTA
153 : integer :: ntemper ! number of temperature points
154 : integer :: use_k_fine ! flag for using fine k-grids for eigenvalues and velocities. 0=no 1=yes
155 : integer :: ep_int_gkk ! flag for interpolate gkk(1) or gamma (0)
156 : integer :: ep_b_min ! first band taken into account in FS integration (if telphint==2)
157 : integer :: ep_b_max ! last band taken into account in FS integration (if telphint==2)
158 : integer :: kptrlatt(3,3) ! kpoint grid generating vectors, as in abinit
159 : integer :: kptrlatt_fine(3,3) ! kpoint grid generating vectors, for fine grid used in FS integration
160 :
161 : real(dp) :: delta_e ! step in electronic energies, around Fermi level
162 : real(dp) :: omega_min,omega_max
163 : real(dp) :: a2fsmear,domega
164 : real(dp) :: nelect ! number of electrons per unit cell, eventually with extra charges for carriers in semiconductors.
165 : real(dp) :: occ_factor ! normalization for integrals over FS, for num of spins, spinors, etc...
166 :
167 : real(dp) :: mustar ! mustar parameter
168 : real(dp) :: fermie ! Fermi energy (Ha), either comes from wfk file or from anaddb input file
169 : real(dp) :: elphsmear ! smearing width for gaussian integration or buffer in energy for
170 : ! calculations with tetrahedra (telphint=0)
171 : real(dp) :: tempermin ! minimum temperature at which resistivity etc are calculated (in K)
172 : real(dp) :: temperinc ! interval temperature grid on which resistivity etc are calculated (in K)
173 :
174 : character(len=fnlen) :: elph_base_name !base name for output files
175 :
176 : integer,allocatable :: qirredtofull(:) !mapping between the qpoints found in the GGK file
177 : !and the array of qpoints generated by the code
178 :
179 : real(dp),allocatable :: wtq(:) !weight for each qpoint in the full grid spqt
180 : !if a point is not in the IBZ ==> wtq=0
181 : !MG we can also use indqpt
182 :
183 : real(dp),allocatable :: n0(:) !DOS at the Fermi level (states/Ha/spin)
184 : real(dp),allocatable :: qpt_full(:,:) !special q points obtained by the Monkhorst & Pack method,
185 : !in reduced coordinates
186 :
187 :
188 : real(dp),allocatable :: gkk_intweight(:,:,:) ! (nFSband,nkpt_fine,nsppol)
189 : !integration weights for gkk matrix elements on FS:
190 : !if ep_keepbands == 0 all are 1
191 : !if ep_keepbands == 1 then = to wtk_phon in elphon
192 : !DOES NOT INCLUDE FACTOR OF 1/nkpt_phon
193 :
194 : real(dp),allocatable :: gkk_velocwtk(:,:,:) ! (nFSband,nkpt_fine,nsppol)
195 :
196 : real(dp),allocatable :: gkk_vvelocwtk(:,:,:) ! (nFSband,nkpt_fine,nsppol)
197 :
198 : real(dp),allocatable :: gkk_qpt(:,:,:,:,:,:) ! (2, ngkkband*ngkkband, nbranch*nbranch, nkpt_phon, nsppol, nqptirred)
199 : !Now gkq contains gkk2 matrices on basic qpts,
200 : !summed over bands if ngkkband==1
201 :
202 :
203 : real(dp),allocatable :: gkk_rpt(:,:,:,:,:,:) ! (2, ngkkband**2, nbranch**2, nkpt_phon, nsppol, nrpt)
204 : !For the moment, gkk_rpt in memory is out of the question
205 : real(dp),allocatable :: gkk2(:,:,:,:,:,:) ! (nbranch, ngkkband,ngkkband, nkpt_phon, nkpt_phon, nsppol)
206 :
207 : real(dp),allocatable :: gamma_qpt(:,:,:,:) !gamma matrices integrated over kpoint coeff
208 : ! and bands: still depends on qpt
209 : ! dims= 2, nbranch**2, nsppol, nqpt
210 : real(dp),allocatable :: gamma_rpt(:,:,:,:)
211 : ! dims= 2, nbranch**2, nsppol, nrpt
212 : !NOTE: choice to put nsppol before or after nqpt is a bit arbitrary
213 : ! abinit uses nband,nkpt,nsppol, but here for convenience nkpt_phon,nsppol,nqpt
214 : ! as interpolation is on qpt
215 :
216 : real(dp),allocatable :: phfrq(:,:) !phonon frequencies
217 : real(dp),allocatable :: a2f(:,:,:) !a2f function
218 :
219 : real(dp),allocatable :: qgrid_data(:,:,:,:) !e-ph values calculated over the irreducible part of the q-grid:
220 : !first entry = index of the q-point,
221 : !second index = branch index
222 : !the third slice contains the frequency, the linewidth and lambda(q,nu)
223 : !for that particular phonon mode
224 : ! dims= nqptirred,elph_ds%nbranch,nsppol,3
225 :
226 : end type elph_type
227 :
228 : public :: elph_ds_clean
229 : !!***
230 :
231 : !----------------------------------------------------------------------
232 :
233 : !!****t* defs_elphon/elph_tr_type
234 : !! NAME
235 : !! elph_tr_type
236 : !!
237 : !! FUNCTION
238 : !! elph_tr_ds contains the necessary data for the transport properties
239 : !!
240 : !! SOURCE
241 :
242 : type,public :: elph_tr_type
243 :
244 : integer :: ifltransport
245 : integer :: unitgkq_trin,unitgkq_trout
246 : integer :: gkqwrite,gkqexist
247 : integer :: onegkksize
248 :
249 : character(len=fnlen) :: ddkfilename
250 :
251 : real(dp),allocatable :: dos_n0(:,:) ! (nT,nsppol) DOS at the Fermi level (states/Ha/spin) at input temperatures
252 : real(dp),allocatable :: dos_n(:,:) ! (nE,nsppol) DOS at the selected energies (states/Ha/spin)
253 : real(dp),allocatable :: en_all(:,:) ! (nE,nsppol) selected energies
254 : real(dp),allocatable :: de_all(:,:) ! (nE,nsppol) differences between selected energies
255 : real(dp),allocatable :: veloc_sq0(:,:,:) ! (3,nsppol,nT)
256 : real(dp),allocatable :: veloc_sq(:,:,:) ! (3,nsppol,nE)
257 :
258 : real(dp),allocatable :: el_veloc(:,:,:,:) ! nkpt nband 3 nsppol
259 : ! the 9 = 3x3 is for the full tensorial transport coefficients
260 : real(dp),allocatable :: gamma_qpt_tr(:,:,:,:,:) ! 2 9 branches**2 nsppol qpt
261 : real(dp),allocatable :: gamma_qpt_trin(:,:,:,:,:) !idem
262 : real(dp),allocatable :: gamma_qpt_trout(:,:,:,:,:) !idem
263 :
264 : real(dp),allocatable :: gamma_rpt_tr(:,:,:,:,:,:,:) !idem
265 : real(dp),allocatable :: gamma_rpt_trin(:,:,:,:,:) !idem
266 : real(dp),allocatable :: gamma_rpt_trout(:,:,:,:,:) !idem
267 :
268 : real(dp),allocatable :: a2f_1d_tr(:,:,:,:,:,:) ! nfreq 9 nsppol 4 n_pair ntemp
269 : real(dp),allocatable :: a2f_1d_trin(:,:,:)
270 : real(dp),allocatable :: a2f_1d_trout(:,:,:)
271 :
272 : real(dp),allocatable :: FSelecveloc_sq(:,:) ! 3 nsppol
273 :
274 : real(dp),allocatable :: tmp_gkk_intweight(:,:,:,:)
275 : real(dp),allocatable :: tmp_gkk_intweight1(:,:,:)
276 : real(dp),allocatable :: tmp_gkk_intweight2(:,:,:)
277 :
278 : real(dp),allocatable :: tmp_velocwtk(:,:,:,:,:)
279 : real(dp),allocatable :: tmp_velocwtk1(:,:,:,:)
280 : real(dp),allocatable :: tmp_velocwtk2(:,:,:,:)
281 :
282 : real(dp),allocatable :: tmp_vvelocwtk(:,:,:,:,:,:)
283 : real(dp),allocatable :: tmp_vvelocwtk1(:,:,:,:,:)
284 : real(dp),allocatable :: tmp_vvelocwtk2(:,:,:,:,:)
285 :
286 : end type elph_tr_type
287 :
288 : public :: elph_tr_ds_clean
289 : !!***
290 :
291 : !----------------------------------------------------------------------
292 :
293 : CONTAINS
294 : !!***
295 :
296 : !----------------------------------------------------------------------
297 :
298 : !!****f* defs_elphon/elph_ds_clean
299 : !!
300 : !! NAME
301 : !! elph_ds_clean
302 : !!
303 : !! FUNCTION
304 : !! deallocate remaining arrays in the elph_ds datastructure
305 : !!
306 : !! INPUTS
307 : !! elph_ds = elphon datastructure
308 : !!
309 : !! SOURCE
310 :
311 15 : subroutine elph_ds_clean(elph_ds)
312 :
313 : !Arguments ------------------------------------
314 : class(elph_type), intent(inout) :: elph_ds
315 : ! *************************************************************************
316 :
317 : !@elph_type
318 15 : ABI_SFREE(elph_ds%qirredtofull)
319 15 : ABI_SFREE(elph_ds%wtq)
320 15 : ABI_SFREE(elph_ds%n0)
321 15 : ABI_SFREE(elph_ds%qpt_full)
322 15 : ABI_SFREE(elph_ds%gkk_intweight)
323 15 : ABI_SFREE(elph_ds%gkk_qpt)
324 15 : ABI_SFREE(elph_ds%gkk_rpt)
325 15 : ABI_SFREE(elph_ds%gkk2)
326 15 : ABI_SFREE(elph_ds%gamma_qpt)
327 15 : ABI_SFREE(elph_ds%gamma_rpt)
328 15 : ABI_SFREE(elph_ds%phfrq)
329 15 : ABI_SFREE(elph_ds%a2f)
330 15 : ABI_SFREE(elph_ds%qgrid_data)
331 :
332 15 : call elph_k_destroy (elph_ds%k_phon)
333 15 : call elph_k_destroy (elph_ds%k_fine)
334 :
335 15 : call elph_ds%k_fine%krank%free()
336 :
337 15 : end subroutine elph_ds_clean
338 : !!***
339 :
340 : !----------------------------------------------------------------------
341 :
342 : !!****f* defs_elphon/elph_tr_ds_clean
343 : !!
344 : !! NAME
345 : !! elph_tr_ds_clean
346 : !!
347 : !! FUNCTION
348 : !! deallocate remaining arrays in the elph_tr_ds datastructure
349 : !!
350 : !! INPUTS
351 : !! elph_tr_ds = elphon transport datastructure
352 : !!
353 : !! SOURCE
354 :
355 15 : subroutine elph_tr_ds_clean(elph_tr_ds)
356 :
357 : !Arguments ------------------------------------
358 : class(elph_tr_type), intent(inout) :: elph_tr_ds
359 : ! *************************************************************************
360 :
361 : !@elph_tr_type
362 15 : ABI_SFREE(elph_tr_ds%el_veloc)
363 15 : ABI_SFREE(elph_tr_ds%FSelecveloc_sq)
364 15 : ABI_SFREE(elph_tr_ds%veloc_sq0)
365 15 : ABI_SFREE(elph_tr_ds%veloc_sq)
366 15 : ABI_SFREE(elph_tr_ds%dos_n0)
367 15 : ABI_SFREE(elph_tr_ds%dos_n)
368 15 : ABI_SFREE(elph_tr_ds%en_all)
369 15 : ABI_SFREE(elph_tr_ds%de_all)
370 15 : ABI_SFREE(elph_tr_ds%gamma_qpt_tr)
371 15 : ABI_SFREE(elph_tr_ds%gamma_qpt_trin)
372 15 : ABI_SFREE(elph_tr_ds%gamma_qpt_trout)
373 15 : ABI_SFREE(elph_tr_ds%gamma_rpt_tr)
374 15 : ABI_SFREE(elph_tr_ds%gamma_rpt_trin)
375 15 : ABI_SFREE(elph_tr_ds%gamma_rpt_trout)
376 15 : ABI_SFREE(elph_tr_ds%a2f_1d_tr)
377 15 : ABI_SFREE(elph_tr_ds%a2f_1d_trin)
378 15 : ABI_SFREE(elph_tr_ds%a2f_1d_trout)
379 15 : ABI_SFREE(elph_tr_ds%tmp_gkk_intweight)
380 15 : ABI_SFREE(elph_tr_ds%tmp_gkk_intweight1)
381 15 : ABI_SFREE(elph_tr_ds%tmp_gkk_intweight2)
382 15 : ABI_SFREE(elph_tr_ds%tmp_velocwtk)
383 15 : ABI_SFREE(elph_tr_ds%tmp_velocwtk1)
384 15 : ABI_SFREE(elph_tr_ds%tmp_velocwtk2)
385 15 : ABI_SFREE(elph_tr_ds%tmp_vvelocwtk)
386 15 : ABI_SFREE(elph_tr_ds%tmp_vvelocwtk1)
387 15 : ABI_SFREE(elph_tr_ds%tmp_vvelocwtk2)
388 :
389 15 : end subroutine elph_tr_ds_clean
390 : !!***
391 :
392 : !----------------------------------------------------------------------
393 :
394 : !!****f* defs_elphon/elph_k_copy
395 : !!
396 : !! NAME
397 : !! elph_k_copy
398 : !!
399 : !! FUNCTION
400 : !! allocate and copy arrays in the elph_k datastructure
401 : !!
402 : !! INPUTS
403 : !! elph_k = elphon k-points datastructure
404 : !!
405 : !! NOTES
406 : !!
407 : !! SOURCE
408 :
409 0 : subroutine elph_k_copy(elph_k_in, elph_k_out)
410 :
411 : !Arguments ------------------------------------
412 : class(elph_kgrid_type), intent(in) :: elph_k_in
413 : class(elph_kgrid_type), intent(out) :: elph_k_out
414 : ! *************************************************************************
415 :
416 : !@elph_kgrid_type
417 0 : elph_k_out%nband = elph_k_in%nband
418 0 : elph_k_out%nsppol = elph_k_in%nsppol
419 0 : elph_k_out%nsym = elph_k_in%nsym
420 :
421 0 : elph_k_out%nkpt = elph_k_in%nkpt
422 0 : elph_k_out%nkptirr = elph_k_in%nkptirr
423 :
424 0 : elph_k_out%my_nkpt = elph_k_in%my_nkpt
425 :
426 0 : ABI_MALLOC(elph_k_out%my_kpt,(elph_k_out%nkpt))
427 0 : elph_k_out%my_kpt = elph_k_in%my_kpt
428 :
429 0 : ABI_MALLOC(elph_k_out%my_ikpt,(elph_k_out%my_nkpt))
430 0 : elph_k_out%my_ikpt = elph_k_in%my_ikpt
431 :
432 0 : ABI_MALLOC(elph_k_out%kptirr,(3,elph_k_out%nkptirr))
433 0 : elph_k_out%kptirr = elph_k_in%kptirr
434 0 : ABI_MALLOC(elph_k_out%wtkirr,(elph_k_out%nkptirr))
435 0 : elph_k_out%wtkirr = elph_k_in%wtkirr
436 :
437 0 : ABI_MALLOC(elph_k_out%wtk,(elph_k_out%nband,elph_k_out%nkpt,elph_k_out%nsppol))
438 0 : elph_k_out%wtk = elph_k_in%wtk
439 0 : ABI_MALLOC(elph_k_out%kpt,(3,elph_k_out%nkpt))
440 0 : elph_k_out%kpt = elph_k_in%kpt
441 :
442 0 : elph_k_out%krank = elph_k_in%krank%copy()
443 :
444 0 : ABI_MALLOC(elph_k_out%irr2full,(elph_k_out%nkptirr))
445 0 : elph_k_out%irr2full = elph_k_in%irr2full
446 0 : ABI_MALLOC(elph_k_out%full2irr,(3,elph_k_out%nkpt))
447 0 : elph_k_out%full2irr = elph_k_in%full2irr
448 0 : ABI_MALLOC(elph_k_out%full2full,(2,elph_k_out%nsym,elph_k_out%nkpt))
449 0 : elph_k_out%full2full = elph_k_in%full2full
450 :
451 0 : ABI_MALLOC(elph_k_out%irredtoGS,(elph_k_out%nkptirr))
452 0 : elph_k_out%irredtoGS = elph_k_in%irredtoGS
453 :
454 0 : end subroutine elph_k_copy
455 : !!***
456 :
457 : !----------------------------------------------------------------------
458 :
459 : !!****f* defs_elphon/elph_k_destroy
460 : !!
461 : !! NAME
462 : !! elph_k_destroy
463 : !!
464 : !! FUNCTION
465 : !! deallocate arrays in the elph_k datastructure
466 : !!
467 : !! INPUTS
468 : !! elph_k = elphon k-points datastructure
469 : !!
470 : !! SOURCE
471 :
472 30 : subroutine elph_k_destroy(elph_k)
473 :
474 : !Arguments ------------------------------------
475 : class(elph_kgrid_type), intent(inout) :: elph_k
476 : ! *************************************************************************
477 :
478 : !@elph_kgrid_type
479 30 : ABI_SFREE(elph_k%irr2full)
480 30 : ABI_SFREE(elph_k%full2irr)
481 30 : ABI_SFREE(elph_k%full2full)
482 30 : ABI_SFREE(elph_k%irredtoGS)
483 30 : ABI_SFREE(elph_k%new_irredtoGS)
484 30 : ABI_SFREE(elph_k%kpt)
485 30 : ABI_SFREE(elph_k%kptirr)
486 30 : ABI_SFREE(elph_k%new_kptirr)
487 30 : ABI_SFREE(elph_k%my_kpt)
488 30 : ABI_SFREE(elph_k%my_ikpt)
489 30 : ABI_SFREE(elph_k%wtk)
490 30 : ABI_SFREE(elph_k%wtq)
491 30 : ABI_SFREE(elph_k%wtkirr)
492 30 : ABI_SFREE(elph_k%new_wtkirr)
493 30 : ABI_SFREE(elph_k%velocwtk)
494 30 : ABI_SFREE(elph_k%vvelocwtk)
495 :
496 30 : call elph_k%krank%free()
497 :
498 30 : end subroutine elph_k_destroy
499 : !!***
500 :
501 : !----------------------------------------------------------------------
502 :
503 : !!****f* defs_elphon/elph_k_procs
504 : !!
505 : !! NAME
506 : !! elph_k_procs
507 : !!
508 : !! FUNCTION
509 : !! allocate kpt to processors, in the elph_k datastructure
510 : !!
511 : !! INPUTS
512 : !! nproc = number of k-parallel processors
513 : !! elph_k = elphon k-points datastructure
514 : !!
515 : !! NOTES
516 : !!
517 : !! SOURCE
518 :
519 15 : subroutine elph_k_procs(nproc, elph_k)
520 :
521 : !Arguments ------------------------------------
522 : !scalars
523 : integer, intent(in) :: nproc
524 : class(elph_kgrid_type), intent(inout) :: elph_k
525 :
526 : integer :: ikpt, me, ik_this_proc
527 : ! *************************************************************************
528 :
529 15 : ABI_SFREE (elph_k%my_kpt)
530 45 : ABI_MALLOC (elph_k%my_kpt, (elph_k%nkpt))
531 :
532 1327 : elph_k%my_kpt = 0
533 15 : elph_k%my_nkpt = 0
534 15 : me = xmpi_comm_rank(xmpi_world)
535 1327 : do ikpt = 1, elph_k%nkpt
536 1312 : elph_k%my_kpt(ikpt) = MOD(ikpt-1, nproc)
537 1327 : if (elph_k%my_kpt(ikpt) == me) elph_k%my_nkpt = elph_k%my_nkpt + 1
538 : end do
539 :
540 : ! create inverse mapping from ik_this_proc to ikpt
541 15 : ABI_SFREE (elph_k%my_ikpt)
542 45 : ABI_MALLOC (elph_k%my_ikpt, (elph_k%my_nkpt))
543 1327 : elph_k%my_ikpt = 0
544 :
545 15 : ik_this_proc = 0
546 1327 : do ikpt = 1, elph_k%nkpt
547 1327 : if (elph_k%my_kpt(ikpt) == me) then
548 1312 : ik_this_proc = ik_this_proc + 1
549 1312 : elph_k%my_ikpt(ik_this_proc) = ikpt
550 : end if
551 : end do
552 15 : ABI_CHECK(ik_this_proc == elph_k%my_nkpt, 'found inconsistent k distribution in processors')
553 :
554 15 : write (std_out,*) 'elph_k_procs : nkpt, distrib = ', elph_k%my_nkpt
555 15 : write (std_out,*) elph_k%my_kpt
556 15 : write (std_out,*) elph_k%my_ikpt
557 :
558 15 : end subroutine elph_k_procs
559 : !!***
560 :
561 : !----------------------------------------------------------------------
562 :
563 : !!****f* defs_elphon/gam_mult_displ
564 : !!
565 : !! NAME
566 : !! gam_mult_displ
567 : !!
568 : !! FUNCTION
569 : !! This routine takes the bare gamma matrices and multiplies them
570 : !! by the displ_red matrices (related to the scalprod variable)
571 : !!
572 : !! INPUTS
573 : !! nbranch = number of phonon branches (3*natom)
574 : !! displ_red = phonon mode displacement vectors in reduced coordinates.
575 : !! gam_bare = bare gamma matrices before multiplication
576 : !!
577 : !! OUTPUT
578 : !! gam_now = output gamma matrices multiplied by displacement matrices
579 : !!
580 : !! SOURCE
581 :
582 48368 : subroutine gam_mult_displ(nbranch, displ_red, gam_bare, gam_now)
583 :
584 : !Arguments -------------------------------
585 : integer, intent(in) :: nbranch
586 : real(dp), intent(in) :: displ_red(2,nbranch,nbranch)
587 : real(dp), intent(in) :: gam_bare(2,nbranch,nbranch)
588 : real(dp), intent(out) :: gam_now(2,nbranch,nbranch)
589 :
590 : !Local variables -------------------------
591 48368 : real(dp) :: zgemm_tmp_mat(2,nbranch,nbranch)
592 : ! *********************************************************************
593 :
594 4337936 : gam_now = zero
595 :
596 48368 : call zgemm('c','n',nbranch,nbranch,nbranch,cone,displ_red,nbranch,gam_bare,nbranch,czero,zgemm_tmp_mat,nbranch)
597 48368 : call zgemm('n','n',nbranch,nbranch,nbranch,cone,zgemm_tmp_mat,nbranch,displ_red,nbranch,czero,gam_now,nbranch)
598 :
599 48368 : end subroutine gam_mult_displ
600 : !!***
601 :
602 : !!****f* ABINIT/complete_gamma
603 : !!
604 : !! NAME
605 : !! complete_gamma
606 : !!
607 : !! FUNCTION
608 : !! Use the set of special q points calculated by the Monkhorst & Pack Technique.
609 : !! Check if all the information for the q points are present in the input gamma matrices.
610 : !! Generate the gamma matrices (already summed over the FS) of the set of q points which
611 : !! samples homogeneously the entire Brillouin zone.
612 : !!
613 : !! INPUTS
614 : !! qpttoqpt = qpoint index mapping under symops
615 : !!
616 : !! OUTPUT
617 : !! gamma_qpt = in/out: set of gamma matrix elements completed and symmetrized
618 : !! gamma_qpt(2,nbranch**2,nsppol,nqpt_full)
619 : !!
620 : !! SOURCE
621 :
622 19 : subroutine complete_gamma(Cryst,nbranch,nsppol,nqptirred,nqpt_full,ep_scalprod,qirredtofull,qpttoqpt,gamma_qpt)
623 :
624 : !Arguments ------------------------------------
625 : !scalars
626 : integer,intent(in) :: nsppol,nbranch,nqptirred,nqpt_full,ep_scalprod
627 : type(crystal_t),intent(in) :: Cryst
628 : !arrays
629 : integer,intent(in) :: qirredtofull(nqptirred)
630 : integer,intent(in) :: qpttoqpt(2,Cryst%nsym,nqpt_full)
631 : real(dp), intent(inout) :: gamma_qpt(2,nbranch**2,nsppol,nqpt_full)
632 :
633 : !Local variables-------------------------------
634 : !scalars
635 : integer :: ibranch,ieqqpt,natom,nsym,iqpt,isppol,isym
636 : integer :: itim,jbranch,neqqpt,iatom,ancestor_iatom,iqpt_fullbz
637 : !arrays
638 38 : integer :: symmetrized_qpt(nqpt_full)
639 38 : integer :: gkk_flag(nbranch,nbranch,nsppol,nqpt_full)
640 : real(dp) :: ss(3,3)
641 38 : real(dp) :: tmp_mat(2,nbranch,nbranch)
642 38 : real(dp) :: tmp_mat2(2,nbranch,nbranch)
643 38 : real(dp) :: ss_allatoms(2,nbranch,nbranch)
644 : complex(dp) :: c_one, c_zero
645 19 : real(dp),allocatable :: gkk_qpt_new(:,:,:),gkk_qpt_tmp(:,:,:)
646 : ! *********************************************************************
647 :
648 19 : c_one = dcmplx(one,zero)
649 19 : c_zero = dcmplx(zero,zero)
650 :
651 19 : natom = Cryst%natom
652 19 : nsym = Cryst%nsym
653 :
654 : !Generation of the gkk matrices relative to the q points
655 : !of the set which samples the entire Brillouin zone
656 :
657 : !set up flags for gamma_qpt matrices we have
658 13163 : gkk_flag = -1
659 125 : do iqpt=1,nqptirred
660 106 : iqpt_fullbz = qirredtofull(iqpt)
661 1902 : gkk_flag(:,:,:,iqpt_fullbz) = 1
662 : end do
663 :
664 899 : symmetrized_qpt(:) = -1
665 :
666 76 : ABI_MALLOC(gkk_qpt_new,(2,nbranch**2,nsppol))
667 57 : ABI_MALLOC(gkk_qpt_tmp,(2,nbranch**2,nsppol))
668 :
669 899 : do iqpt=1,nqpt_full
670 : !
671 : ! Already symmetrized?
672 880 : if (symmetrized_qpt(iqpt) == 1) cycle
673 :
674 4130 : gkk_qpt_new(:,:,:) = zero
675 :
676 : ! loop over qpoints equivalent to iqpt
677 : neqqpt=0
678 : ! do not use time reversal symmetry to complete the qpoints:
679 : ! do not know what happens to the gamma matrices
680 : ! 11/2011: MJV: time reversal is needed here if inversion is absent
681 : ! - used in read_gkk and all reductions of q-points by symmetry.
682 :
683 318 : do itim=1,2
684 9822 : do isym=1,nsym
685 : ! ieqqpt is sent onto iqpt by itim/isym
686 9504 : ieqqpt = qpttoqpt(itim,isym,iqpt)
687 :
688 :
689 9504 : if (gkk_flag(1,1,1,ieqqpt) == -1) cycle
690 : ! if we have information on this qpt
691 : ! iqpt is equivalent to ieqqpt: get it from file or memory
692 121932 : gkk_qpt_tmp(:,:,:) = gamma_qpt(:,:,:,ieqqpt)
693 :
694 3164 : neqqpt=neqqpt+1
695 :
696 : !
697 : ! MJV note 02/2010:
698 : ! the correspondence of symrel and symrec in the different cases, symmetrizing there
699 : ! and back, has been fixed in the cases with and without scalprod (ie cartesian
700 : ! and reduced real space coordinates) with respect to a calculation with no symmetries
701 : ! I believe everything is settled, but still do not know why the 2 versions of the ss
702 : ! matrices here use different rel/rec, instead of just being multiplied by the rprim gprim...
703 : !
704 3164 : if (ep_scalprod==1) then
705 : ! ss(ii,jj)=ss(ii,jj)+Cryst%rprimd(ii,kk)*Cryst%symrel(kk,ll,isym)*Cryst%gprimd(ll,jj)
706 63248 : ss(:,:) = MATMUL(Cryst%rprimd, MATMUL(Cryst%symrel(:,:,isym), Cryst%gprimd))
707 : else
708 : ! ss(ii,jj) = Cryst%symrec(ii,jj,isym)
709 34996 : ss(:,:) = Cryst%symrec(:,:,isym)
710 : end if
711 :
712 124964 : ss_allatoms(:,:,:) = zero
713 6648 : do iatom=1,natom
714 3484 : ancestor_iatom = Cryst%indsym(4,isym,iatom)
715 48456 : ss_allatoms(1, (ancestor_iatom-1)*3+1:(ancestor_iatom-1)*3+3, (iatom-1)*3+1:(iatom-1)*3+3) = ss(1:3,1:3)
716 : end do
717 :
718 :
719 : ! NOTE ssinv(ii,jj)=ssinv(ii,jj)+Cryst%gprimd(ii,kk)*rprimd(jj,ll)*Cryst%symrec(ll,kk,isym)
720 :
721 6692 : do isppol=1,nsppol
722 : ! multiply by the ss matrices
723 129676 : tmp_mat2(:,:,:) = zero
724 13264 : tmp_mat(:,:,:) = reshape(gkk_qpt_tmp(:,:,isppol),(/2,nbranch,nbranch/))
725 :
726 : call ZGEMM ('N','N',nbranch,nbranch,nbranch,&
727 3316 : & c_one,ss_allatoms,nbranch,tmp_mat,nbranch,c_zero,tmp_mat2,nbranch)
728 :
729 : call ZGEMM ('N','T',nbranch,nbranch,nbranch,&
730 3316 : & c_one,tmp_mat2,nbranch,ss_allatoms,nbranch,c_zero,tmp_mat,nbranch)
731 :
732 : ! add to gkk_qpt_new
733 23728 : do ibranch =1,nbranch
734 52708 : do jbranch =1,nbranch
735 : gkk_qpt_new(:,(jbranch-1)*nbranch+ibranch,isppol) = &
736 126360 : & gkk_qpt_new(:,(jbranch-1)*nbranch+ibranch,isppol) + tmp_mat(:,jbranch,ibranch)
737 : end do
738 : end do
739 : end do ! isppol
740 : !
741 : end do ! isym
742 : end do ! itim
743 : !
744 106 : ABI_CHECK(neqqpt>0,'no q-points found equivalent to iqpt ')
745 : ! Divide by number of equivalent qpts found.
746 4130 : gkk_qpt_new(:,:,:) = gkk_qpt_new(:,:,:)/neqqpt
747 :
748 : ! copy the symmetrized version into all the equivalent qpoints, appropriately transformed
749 337 : do itim=1,2
750 10596 : do isym=1,nsym
751 : ! ieqqpt is sent onto iqpt by itim/isym
752 9504 : ieqqpt = qpttoqpt(itim,isym,iqpt)
753 :
754 9504 : if (symmetrized_qpt(ieqqpt) /= -1) cycle
755 27688 : gkk_qpt_tmp(:,:,:) = zero
756 :
757 : ! use symrec matrices to get inverse transform from isym^{-1}
758 880 : if (ep_scalprod==1) then
759 : ! Use inverse of symop matrix here to get back to ieqqpt (inv+transpose is in symrec and in gprimd)
760 : ! ss(ii,jj)=ss(ii,jj)+Cryst%rprimd(ii,kk)*Cryst%symrec(ll,kk,isym)*Cryst%gprimd(ll,jj)
761 4288 : ss(:,:) = MATMUL(Cryst%rprimd, MATMUL(TRANSPOSE(Cryst%symrec(:,:,isym)), Cryst%gprimd))
762 : else
763 : ! ss(ii,jj) = Cryst%symrel(jj,ii,isym)
764 11024 : ss(:,:) = TRANSPOSE(Cryst%symrel(:,:,isym))
765 : end if
766 :
767 29296 : ss_allatoms(:,:,:) = zero
768 1784 : do iatom=1,natom
769 904 : ancestor_iatom = Cryst%indsym(4,isym,iatom)
770 12632 : ss_allatoms(1, (ancestor_iatom-1)*3+1:(ancestor_iatom-1)*3+3, (iatom-1)*3+1:(iatom-1)*3+3) = ss(1:3,1:3)
771 : end do
772 :
773 : ! ! Use inverse of symop matrix here to get back to ieqqpt
774 : ! ssinv(ii,jj)=ssinv(ii,jj)+gprimd(ii,kk)*rprimd(jj,ll)*Cryst%symrel(kk,ll,isym)
775 :
776 1768 : do isppol=1,nsppol
777 : ! multiply by the ss^{-1} matrices
778 29544 : tmp_mat2(:,:,:) = zero
779 3552 : tmp_mat(:,:,:) = reshape(gkk_qpt_new(:,:,isppol),(/2,nbranch,nbranch/))
780 :
781 : call ZGEMM ('N','N',nbranch,nbranch,nbranch,&
782 888 : & c_one,ss_allatoms,nbranch,tmp_mat,nbranch,c_zero,tmp_mat2,nbranch)
783 :
784 : call ZGEMM ('N','T',nbranch,nbranch,nbranch,&
785 888 : & c_one,tmp_mat2,nbranch,ss_allatoms,nbranch,c_zero,tmp_mat,nbranch)
786 :
787 : ! FIXME: the following could just be a reshape
788 3624 : do ibranch =1,nbranch
789 12264 : do jbranch =1,nbranch
790 : gkk_qpt_tmp(:,(jbranch-1)*nbranch+ibranch,isppol) =&
791 28656 : & tmp_mat(:,jbranch,ibranch)
792 : end do
793 : end do
794 11476 : if (gkk_flag (1,1,isppol,ieqqpt) == -1) gkk_flag (:,:,isppol,ieqqpt) = 0
795 : end do ! end isppol do
796 :
797 : ! save symmetrized matrices for qpt ieqqpt
798 27688 : gamma_qpt(:,:,:,ieqqpt) = gkk_qpt_tmp(:,:,:)
799 :
800 9716 : symmetrized_qpt(ieqqpt) = 1
801 :
802 : end do !isym
803 : end do !itim
804 : end do !iqpt
805 :
806 19 : ABI_FREE(gkk_qpt_new)
807 19 : ABI_FREE(gkk_qpt_tmp)
808 :
809 19 : end subroutine complete_gamma
810 : !!***
811 :
812 : !!****f* ABINIT/complete_gamma_tr
813 : !!
814 : !! NAME
815 : !! complete_gamma_tr
816 : !!
817 : !! FUNCTION
818 : !! Use the set of special q points calculated by the Monkhorst & Pack Technique.
819 : !! Check if all the information for the q points are present in
820 : !! the input gamma transport matrices.
821 : !! Generate the gamma transport matrices (already summed over the FS) of the set of q points which
822 : !! samples homogeneously the entire Brillouin zone.
823 : !!
824 : !! INPUTS
825 : !! crystal<crystal_t>=data type gathering info on the crystalline structure.
826 : !! ep_scalprod= flag for scalar product of gkk with phonon displacement vectors
827 : !! nbranch=number of phonon branches = 3*natom
828 : !! nqptirred=nqpt irred BZ
829 : !! nqpt_full=nqpt full BZ
830 : !! nsppol=number of spins
831 : !! qirredtofull= mapping irred to full qpoints
832 : !! qpttoqpt = qpoint index mapping under symops
833 : !!
834 : !! OUTPUT
835 : !! gamma_qpt_tr = in/out: set of gamma matrix elements completed and symmetrized
836 : !! gamma_qpt_tr(2,9,nbranch*nbranch,nsppol,nqpt_full)
837 : !!
838 : !! SOURCE
839 :
840 10 : subroutine complete_gamma_tr(crystal,ep_scalprod,nbranch,nqptirred,nqpt_full,nsppol,gamma_qpt_tr,qirredtofull,qpttoqpt)
841 :
842 : use m_linalg_interfaces
843 :
844 : !Arguments ------------------------------------
845 : !scalars
846 : integer, intent(in) :: nbranch,nqptirred,nqpt_full,nsppol, ep_scalprod
847 : type(crystal_t),intent(in) :: crystal
848 : !arrays
849 : integer,intent(in) :: qpttoqpt(2,crystal%nsym,nqpt_full)
850 : integer,intent(in) :: qirredtofull(nqptirred)
851 : real(dp), intent(inout) :: gamma_qpt_tr(2,9,nbranch*nbranch,nsppol,nqpt_full)
852 :
853 : !Local variables-------------------------------
854 : !scalars
855 : integer :: ieqqpt,iqpt,isppol,isym
856 : integer :: itim,neqqpt
857 : integer :: iatom,ancestor_iatom
858 : integer :: iqpt_fullbz,imode, itensor,reim
859 : real(dp),parameter :: tol=2.d-8
860 : !arrays
861 20 : integer :: symrel(3,3,crystal%nsym),symrec(3,3,crystal%nsym)
862 20 : integer :: symmetrized_qpt(nqpt_full)
863 20 : integer :: gkk_flag(nbranch,nbranch,nsppol,nqpt_full)
864 : real(dp) :: gprimd(3,3),rprimd(3,3)
865 : real(dp) :: ss(3,3), sscart(3,3)
866 20 : real(dp) :: tmp_mat(nbranch,nbranch)
867 20 : real(dp) :: tmp_mat2(nbranch,nbranch)
868 : real(dp) :: tmp_tensor(3,3)
869 : real(dp) :: tmp_tensor2(3,3)
870 20 : real(dp) :: ss_allatoms(nbranch,nbranch)
871 10 : real(dp),allocatable :: gkk_qpt_new(:,:,:,:),gkk_qpt_tmp(:,:,:,:)
872 : ! *********************************************************************
873 :
874 130 : gprimd = crystal%gprimd
875 130 : rprimd = crystal%rprimd
876 :
877 4066 : symrec = crystal%symrec
878 4066 : symrel = crystal%symrel
879 :
880 : !Generation of the gkk matrices relative to the q points
881 : !of the set which samples the entire Brillouin zone
882 :
883 : !set up flags for gamma_qpt matrices we have
884 2570 : gkk_flag = -1
885 46 : do iqpt=1,nqptirred
886 36 : iqpt_fullbz = qirredtofull(iqpt)
887 1234 : gkk_flag(:,:,:,iqpt_fullbz) = 1
888 : end do
889 :
890 90 : symmetrized_qpt(:) = -1
891 : ! isppol=1
892 :
893 40 : ABI_MALLOC(gkk_qpt_new,(2,9,nbranch*nbranch, nsppol))
894 30 : ABI_MALLOC(gkk_qpt_tmp,(2,9,nbranch*nbranch, nsppol))
895 :
896 90 : do iqpt=1,nqpt_full
897 :
898 : ! Already symmetrized?
899 80 : if (symmetrized_qpt(iqpt) == 1) cycle
900 :
901 27288 : gkk_qpt_new(:,:,:,:) = zero
902 :
903 : ! loop over qpoints equivalent to iqpt
904 : neqqpt=0
905 : ! do not use time reversal symmetry to complete the qpoints:
906 : ! do not know what happens to the gamma matrices
907 :
908 108 : do itim=1,2
909 2220 : do isym=1,crystal%nsym
910 : ! ieqqpt is sent onto iqpt by itim/isym
911 2112 : ieqqpt = qpttoqpt(itim,isym,iqpt)
912 :
913 2112 : if (gkk_flag(1,1,1,ieqqpt) == -1) cycle
914 : ! if we have information on this qpt
915 : ! iqpt is equivalent to ieqqpt: get it from file or memory
916 800832 : gkk_qpt_tmp(:,:,:,:) = gamma_qpt_tr(:,:,:,:,ieqqpt)
917 :
918 1248 : neqqpt=neqqpt+1
919 :
920 : !
921 : ! MJV note 02/2010:
922 : ! the correspondence of symrel and symrec in the different cases, symmetrizing there
923 : ! and back, has been fixed in the cases with and without scalprod (ie cartesian
924 : ! and reduced real space coordinates) with respect to a calculation with no symmetries
925 : ! I believe everything is settled, but still do not know why the 2 versions of the ss
926 : ! matrices here use different rel/rec, instead of just being multiplied by the rprim gprim...
927 : !
928 : ! sscart(ii,jj)=sscart(ii,jj)+rprimd(ii,kk)*symrel(kk,ll,isym)*gprimd(ll,jj)
929 113568 : sscart(:,:) = MATMUL(rprimd, MATMUL(symrel(:,:,isym), gprimd))
930 1248 : if (ep_scalprod==1) then
931 : ! ss(ii,jj)=ss(ii,jj)+rprimd(ii,kk)*symrel(kk,ll,isym)*gprimd(ll,jj)
932 85904 : ss(:,:) = MATMUL(rprimd, MATMUL(symrel(:,:,isym), gprimd))
933 : else
934 : ! ss(ii,jj) = symrec(ii,jj,isym)
935 3952 : ss(:,:) = symrec(:,:,isym)
936 : end if
937 :
938 35424 : ss_allatoms(:,:) = zero
939 3136 : do iatom=1,crystal%natom
940 1888 : ancestor_iatom = crystal%indsym(4,isym,iatom)
941 : ss_allatoms((ancestor_iatom-1)*3+1:(ancestor_iatom-1)*3+3,&
942 25792 : & (iatom-1)*3+1: (iatom-1)*3+3) = ss(1:3,1:3)
943 : end do
944 :
945 :
946 : ! NOTE ssinv(ii,jj)=ssinv(ii,jj)+gprimd(ii,kk)*rprimd(jj,ll)*symrec(ll,kk,isym)
947 :
948 2568 : do isppol=1,nsppol
949 :
950 : ! for each tensor component, rotate the cartesian directions of phonon modes
951 12480 : do itensor = 1, 9
952 34944 : do reim=1,2 ! Real/Imaginary parts
953 : ! Multiply by the ss matrices
954 637632 : tmp_mat2(:,:) = zero
955 67392 : tmp_mat(:,:) = reshape(gkk_qpt_tmp(reim,itensor,:,isppol),(/nbranch,nbranch/))
956 22464 : call DGEMM ('N','N',nbranch,nbranch,nbranch,one,ss_allatoms,nbranch,tmp_mat,nbranch,zero,tmp_mat2,nbranch)
957 22464 : call DGEMM ('N','T',nbranch,nbranch,nbranch,one,tmp_mat2,nbranch,ss_allatoms,nbranch,zero,tmp_mat,nbranch)
958 569376 : gkk_qpt_tmp(reim,itensor,:,isppol) = reshape (tmp_mat, (/nbranch*nbranch/))
959 : enddo
960 : end do ! itensor
961 :
962 : ! for each cartesian direction/phonon mode, rotate the tensor components
963 29760 : do imode = 1, nbranch*nbranch
964 86784 : do reim=1,2 ! Real/Imaginary parts
965 57024 : tmp_tensor2(:,:) = zero
966 57024 : tmp_tensor(:,:) = reshape(gkk_qpt_tmp(reim,:,imode,isppol),(/3,3/))
967 57024 : call DGEMM ('N','N',3,3,3,one,sscart,3,tmp_tensor,3,zero,tmp_tensor2,3)
968 57024 : call DGEMM ('N','T',3,3,3,one,tmp_tensor2,3,sscart,3,zero,tmp_tensor,3)
969 598752 : gkk_qpt_tmp(reim,:,imode,isppol) = reshape (tmp_tensor, (/9/)) ! modified by BX
970 : enddo
971 : end do ! imode
972 :
973 : ! add to gkk_qpt_new
974 800832 : gkk_qpt_new(:,:,:,isppol) = gkk_qpt_new(:,:,:,isppol) + gkk_qpt_tmp(:,:,:,isppol)
975 :
976 : end do ! end isppol do
977 :
978 : end do ! end isym do
979 : end do ! end itim do
980 :
981 36 : ABI_CHECK(neqqpt>0,'no q-points found equivalent to iqpt ')
982 :
983 : ! divide by number of equivalent qpts found
984 27288 : gkk_qpt_new = gkk_qpt_new/neqqpt
985 :
986 :
987 : ! copy the symmetrized version into all the equivalent qpoints, appropriately transformed
988 118 : do itim=1,2
989 2264 : do isym=1,crystal%nsym
990 : ! ieqqpt is sent onto iqpt by itim/isym
991 2112 : ieqqpt = qpttoqpt(itim,isym,iqpt)
992 :
993 2112 : if (symmetrized_qpt(ieqqpt) /= -1) cycle
994 56608 : gkk_qpt_tmp = zero
995 :
996 : ! use symrec matrices to get inverse transform from isym^{-1}
997 : ! Use inverse of symop matrix here to get back to ieqqpt (inv+transpose is in symrec and in gprimd)
998 : ! sscart(ii,jj)=sscart(ii,jj)+rprimd(ii,kk)*symrec(ll,kk,isym)*gprimd(ll,jj)
999 : ! sscart(ii,jj)=sscart(ii,jj)+rprimd(ii,kk)*symrec(ll,kk,isym)*gprimd(ll,jj)
1000 7280 : sscart(:,:) = MATMUL(rprimd, MATMUL(TRANSPOSE(symrec(:,:,isym)), gprimd))
1001 80 : if (ep_scalprod==1) then
1002 : ! Use inverse of symop matrix here to get back to ieqqpt (inv+transpose is in symrec and in gprimd)
1003 : ! ss(ii,jj)=ss(ii,jj)+rprimd(ii,kk)*symrec(ll,kk,isym)*gprimd(ll,jj)
1004 5824 : ss(:,:) = MATMUL(rprimd, MATMUL(TRANSPOSE(symrec(:,:,isym)), gprimd))
1005 :
1006 : else
1007 : ! ss(ii,jj) = symrel(jj,ii,isym)
1008 208 : ss(:,:) = TRANSPOSE(symrel(:,:,isym))
1009 : end if
1010 :
1011 2480 : ss_allatoms(:,:) = zero
1012 208 : do iatom=1,crystal%natom
1013 128 : ancestor_iatom = crystal%indsym(4,isym,iatom)
1014 : ss_allatoms((ancestor_iatom-1)*3+1:(ancestor_iatom-1)*3+3,&
1015 1744 : & (iatom-1)*3+1: (iatom-1)*3+3) = ss(1:3,1:3)
1016 : end do
1017 :
1018 : ! ! Use inverse of symop matrix here to get back to ieqqpt
1019 : ! ssinv(ii,jj)=ssinv(ii,jj)+gprimd(ii,kk)*rprimd(jj,ll)*symrel(kk,ll,isym)
1020 :
1021 160 : do isppol=1,nsppol
1022 800 : do itensor = 1, 9
1023 2240 : do reim=1,2 ! Real/Imaginary parts
1024 : ! Multiply by the ss^{-1} matrices
1025 44640 : tmp_mat2(:,:) = zero
1026 4320 : tmp_mat(:,:) = reshape(gkk_qpt_new(reim,itensor,:,isppol),(/nbranch,nbranch/))
1027 1440 : call DGEMM ('N','N',nbranch,nbranch,nbranch,one,ss_allatoms,nbranch,tmp_mat,nbranch,zero,tmp_mat2,nbranch)
1028 1440 : call DGEMM ('N','T',nbranch,nbranch,nbranch,one,tmp_mat2,nbranch,ss_allatoms,nbranch,zero,tmp_mat,nbranch)
1029 39888 : gkk_qpt_tmp(reim,itensor,:,isppol) = reshape (tmp_mat, (/nbranch*nbranch/))
1030 : enddo
1031 : end do ! itensor
1032 :
1033 : ! for each cartesian direction/phonon mode, rotate the tensor components
1034 2096 : do imode = 1, nbranch*nbranch
1035 6128 : do reim=1,2 ! Real/Imaginary parts
1036 4032 : tmp_tensor2(:,:) = zero
1037 4032 : tmp_tensor(:,:) = reshape(gkk_qpt_tmp(reim,:,imode,isppol),(/3,3/))
1038 4032 : call DGEMM ('N','N',3,3,3,one,sscart,3,tmp_tensor,3,zero,tmp_tensor2,3)
1039 4032 : call DGEMM ('N','T',3,3,3,one,tmp_tensor2,3,sscart,3,zero,tmp_tensor,3)
1040 : ! gkk_qpt_new(:,:,imode,isppol) = reshape (tmp_tensor, (/2,9/)) ! Modified by BX
1041 42336 : gkk_qpt_tmp(reim,:,imode,isppol) = reshape (tmp_tensor, (/9/)) ! Modified by BX
1042 : enddo
1043 : end do ! imode
1044 :
1045 160 : if (gkk_flag (1,1,isppol,ieqqpt) == -1) then
1046 1292 : gkk_flag (:,:,isppol,ieqqpt) = 0
1047 : end if
1048 :
1049 : end do ! end isppol do
1050 :
1051 :
1052 : ! save symmetrized matrices for qpt ieqqpt
1053 56608 : gamma_qpt_tr(:,:,:,:,ieqqpt) = gkk_qpt_tmp(:,:,:,:)
1054 :
1055 152 : symmetrized_qpt(ieqqpt) = 1
1056 :
1057 : end do ! end isym do
1058 : end do ! end itim do
1059 :
1060 : end do
1061 : !end iqpt do
1062 :
1063 10 : ABI_FREE(gkk_qpt_new)
1064 10 : ABI_FREE(gkk_qpt_tmp)
1065 :
1066 10 : end subroutine complete_gamma_tr
1067 : !!***
1068 :
1069 : !----------------------------------------------------------------------
1070 :
1071 : !!****f* m_fstab/mkqptequiv
1072 : !! NAME
1073 : !! mkqptequiv
1074 : !!
1075 : !! FUNCTION
1076 : !! This routine determines the equivalence between
1077 : !! 1) qpoints and fermi surface kpoints
1078 : !! 2) qpoints under symmetry operations
1079 : !!
1080 : !! INPUTS
1081 : !! Cryst<crystal_t>=Info on unit cell and symmetries.
1082 : !! kpt_phon = fermi surface kpoints
1083 : !! nkpt_phon = number of kpoints in the full FS set
1084 : !! nqpt = number of qpoints
1085 : !! qpt_full = qpoint coordinates
1086 : !!
1087 : !! OUTPUT
1088 : !! FSfullpqtofull = mapping of k + q onto k' for k and k' in full BZ
1089 : !! qpttoqpt(itim,isym,iqpt) = qpoint index which transforms to iqpt under isym and with time reversal itim.
1090 : !!
1091 : !! NOTES
1092 : !! REMOVED 3/6/2008: much too large matrix, and not used at present
1093 : !! FStoqpt = mapping of kpoint pairs (1 irreducible and 1 full) to qpoints
1094 : !!
1095 : !! SOURCE
1096 :
1097 15 : subroutine mkqptequiv(FSfullpqtofull,Cryst,kpt_phon,nkpt_phon,nqpt,qpttoqpt,qpt_full,mqtofull)
1098 :
1099 : !Arguments ------------------------------------
1100 : !scalars
1101 : integer,intent(in) :: nkpt_phon,nqpt
1102 : type(crystal_t),intent(in) :: Cryst
1103 : !arrays
1104 : integer,intent(out) :: FSfullpqtofull(nkpt_phon,nqpt),qpttoqpt(2,Cryst%nsym,nqpt)
1105 : integer,intent(out),optional :: mqtofull(nqpt)
1106 : real(dp),intent(in) :: kpt_phon(3,nkpt_phon),qpt_full(3,nqpt)
1107 :
1108 : !Local variables-------------------------------
1109 : !scalars
1110 : integer :: ikpt_phon,iFSqpt,iqpt,isym,symrankkpt_phon
1111 : !character(len=500) :: message
1112 15 : type(krank_t) :: krank
1113 : !arrays
1114 : real(dp) :: tmpkpt(3),gamma_kpt(3)
1115 : ! *************************************************************************
1116 :
1117 15 : call wrtout(std_out,' mkqptequiv : making rankkpt_phon and invrankkpt_phon',"COLL")
1118 :
1119 15 : call krank%init(nkpt_phon, kpt_phon)
1120 :
1121 10631 : FSfullpqtofull = -999
1122 : gamma_kpt(:) = zero
1123 :
1124 1327 : do ikpt_phon=1,nkpt_phon
1125 11823 : do iqpt=1,nqpt
1126 : ! tmpkpt = jkpt = ikpt + qpt
1127 41984 : tmpkpt(:) = kpt_phon(:,ikpt_phon) + qpt_full(:,iqpt)
1128 :
1129 : ! which kpt is it among the full FS kpts?
1130 10496 : symrankkpt_phon = krank%get_rank(tmpkpt)
1131 :
1132 10496 : FSfullpqtofull(ikpt_phon,iqpt) = krank%invrank(symrankkpt_phon)
1133 11808 : if (FSfullpqtofull(ikpt_phon, iqpt) == -1) then
1134 0 : ABI_ERROR("looks like no kpoint equiv to k+q !!!")
1135 : end if
1136 :
1137 : end do
1138 : end do
1139 :
1140 15 : if (present(mqtofull)) then
1141 0 : do iqpt=1,nqpt
1142 0 : tmpkpt(:) = gamma_kpt(:) - qpt_full(:,iqpt)
1143 :
1144 : ! which kpt is it among the full FS kpts?
1145 0 : symrankkpt_phon = krank%get_rank(tmpkpt)
1146 :
1147 0 : mqtofull(iqpt) = krank%invrank(symrankkpt_phon)
1148 0 : if (mqtofull(iqpt) == -1) then
1149 0 : ABI_ERROR("looks like no kpoint equiv to -q !!!")
1150 : end if
1151 : end do
1152 : end if
1153 :
1154 15 : call krank%free()
1155 :
1156 : ! start over with q grid
1157 15 : call wrtout(std_out,' mkqptequiv : FSfullpqtofull made. Do qpttoqpt',"COLL")
1158 :
1159 15 : call krank%init(nqpt, qpt_full)
1160 :
1161 15399 : qpttoqpt(:,:,:) = -1
1162 135 : do iFSqpt=1,nqpt
1163 5223 : do isym=1,Cryst%nsym
1164 : tmpkpt(:) = Cryst%symrec(:,1,isym)*qpt_full(1,iFSqpt) &
1165 : + Cryst%symrec(:,2,isym)*qpt_full(2,iFSqpt) &
1166 20352 : + Cryst%symrec(:,3,isym)*qpt_full(3,iFSqpt)
1167 :
1168 5088 : symrankkpt_phon = krank%get_rank(tmpkpt)
1169 5088 : if (krank%invrank(symrankkpt_phon) == -1) then
1170 0 : ABI_ERROR("looks like no kpoint equiv to q by symmetry without time reversal!!!")
1171 : end if
1172 5088 : qpttoqpt(1,isym,krank%invrank(symrankkpt_phon)) = iFSqpt
1173 :
1174 20352 : tmpkpt = -tmpkpt
1175 5088 : symrankkpt_phon = krank%get_rank(tmpkpt)
1176 5088 : if (krank%invrank(symrankkpt_phon) == -1) then
1177 0 : ABI_ERROR('looks like no kpoint equiv to q by symmetry with time reversal!!!')
1178 : end if
1179 5208 : qpttoqpt(2,isym,krank%invrank(symrankkpt_phon)) = iFSqpt
1180 : end do
1181 : end do
1182 :
1183 15 : call krank%free()
1184 :
1185 15 : end subroutine mkqptequiv
1186 : !!***
1187 :
1188 504 : end module defs_elphon
1189 : !!***
|
