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