Line data Source code
1 : !!****m* ABINIT/m_eig2d
2 : !! NAME
3 : !! m_eig2d
4 : !!
5 : !! FUNCTION
6 : !! This module contains utilities to analyze and retrieve information
7 : !! from the second order derivative of the eigen-energies wrt
8 : !! displacements.
9 : !!
10 : !! COPYRIGHT
11 : !! Copyright (C) 2014-2026 ABINIT group (SP, PB, XG)
12 : !! This file is distributed under the terms of the
13 : !! GNU General Public License, see ~abinit/COPYING
14 : !! or http://www.gnu.org/copyleft/gpl.txt .
15 : !!
16 : !! SOURCE
17 :
18 : #if defined HAVE_CONFIG_H
19 : #include "config.h"
20 : #endif
21 :
22 : #include "abi_common.h"
23 :
24 : MODULE m_eig2d
25 :
26 : use defs_basis
27 : use m_errors
28 : use m_abicore
29 : use m_nctk
30 : use netcdf
31 : use m_xmpi
32 : !use m_ebands
33 : use m_cgtools
34 : use m_hdr
35 : use m_dtset
36 : use m_dtfil
37 : use m_ddb_hdr
38 : use m_ddb
39 :
40 : use defs_datatypes, only : pseudopotential_type
41 : use defs_abitypes, only : MPI_type
42 : use m_time, only : timab
43 : use m_fstrings, only : strcat
44 : use m_crystal, only : crystal_t
45 : use m_pawtab, only : pawtab_type
46 : use m_double_grid,only : kptfine_av
47 : use m_mpinfo, only : distrb2, proc_distrb_cycle
48 : use m_ebands, only : ebands_t
49 :
50 : implicit none
51 :
52 : private
53 :
54 : public :: eigr2d_init ! Main creation method of EIG2D.nc files.
55 : public :: eigr2d_ncwrite ! Dump the object into NETCDF file.
56 : public :: eigr2d_free ! Destruction method.
57 : public :: fan_init ! Main creation method of Fan.nc files.
58 : public :: fan_ncwrite ! Dump the object into NETCDF file.
59 : public :: fan_free ! Destruction method.
60 : public :: gkk_init ! Main creation method of GKK.nc files.
61 : public :: gkk_ncwrite ! Dump the object into NETCDF file.
62 : public :: gkk_free ! Destruction method.
63 :
64 : public :: eig2tot ! This routine calculates the second-order eigenvalues.
65 : public :: outbsd ! output bsd file for one perturbation (used for elphon calculations in anaddb)
66 : public :: eig2stern
67 : public :: elph2_fanddw ! Calculates the zero-point motion corrections
68 :
69 : !!***
70 :
71 : !!****t* m_eig2d/eigr2d_t
72 : !! NAME
73 : !! eig2d_t
74 : !!
75 : !! FUNCTION
76 : !! It contains information about the second-order derivative of the
77 : !! eigenenergies wrt atomic displacement
78 : !!
79 : !! SOURCE
80 :
81 : type,public :: eigr2d_t
82 :
83 : ! WARNING : if you modify this datatype, please check whether there might be
84 : ! creation/destruction/copy routines,
85 : ! declared in another part of ABINIT, that might need to take into account your
86 : ! modification.
87 :
88 : integer :: mband ! Max number of bands i.e MAXVAL(nband) (to dimension arrays)
89 : integer :: nsppol ! number of spin-polarization
90 : integer :: nkpt ! number of k points
91 : integer :: natom ! number of atoms
92 :
93 : real(dp),allocatable :: eigr2d(:,:,:,:,:,:,:)
94 : ! eigr2d(2,mband*nsppol,nkpt,3,natom,3,natom)
95 : ! Second-order derivative of eigenergies (real,im) at each
96 : ! spin,band,k-point,dir1,dir2,natom1,natom2 .
97 :
98 :
99 : end type eigr2d_t
100 : !!***
101 :
102 : !!****t* m_eig2d/fan_t
103 : !! NAME
104 : !! fan_t
105 : !!
106 : !! FUNCTION
107 : !! It contains information about the second-order derivative of the
108 : !! eigenenergies wrt atomic displacement
109 : !!
110 : !! SOURCE
111 :
112 : type,public :: fan_t
113 :
114 : ! WARNING : if you modify this datatype, please check whether there might be
115 : ! creation/destruction/copy routines,
116 : ! declared in another part of ABINIT, that might need to take into account your
117 : ! modification.
118 :
119 : integer :: mband ! Max number of bands i.e MAXVAL(nband) (to dimension arrays)
120 : integer :: nsppol ! number of spin-polarization
121 : integer :: nkpt ! number of k points
122 : integer :: natom ! number of atoms
123 :
124 : real(dp),allocatable :: fan2d(:,:,:,:,:,:,:)
125 : ! fan2d(2*mband*nsppol,nkpt,3,natom,3,natom,mband)
126 : ! Second-order derivative of the eigenergies (real,im) at each
127 : ! ispin,iband(real,im),k-point,dir1,dir2,natom1,natom2,jband
128 :
129 : end type fan_t
130 : !!***
131 :
132 : !!****t* m_eig2d/gkk_t
133 : !! NAME
134 : !! gkk_t
135 : !!
136 : !! FUNCTION
137 : !! It contains information about the second-order derivative of the
138 : !! eigenenergies wrt atomic displacement
139 : !!
140 : !! SOURCE
141 :
142 : type,public :: gkk_t
143 :
144 : ! WARNING : if you modify this datatype, please check whether there might be
145 : ! creation/destruction/copy routines,
146 : ! declared in another part of ABINIT, that might need to take into account your
147 : ! modification.
148 :
149 : integer :: mband ! Max number of bands i.e MAXVAL(nband) (to dimension arrays)
150 : integer :: nsppol ! number of spin-polarization
151 : integer :: nkpt ! number of k points
152 : integer :: natom ! number of atoms
153 : integer :: ncart ! number of cartesian directions
154 :
155 : real(dp),allocatable :: gkk2d(:,:,:,:,:)
156 : ! gkk2d(2*mband*nsppol,nkpt,ncart,natom,mband)
157 : ! Second-order derivative of the eigenergies (real,im) at each
158 : ! ispin,iband(real,im),k-point,dir1,natom1,jband
159 :
160 : end type gkk_t
161 : !!***
162 :
163 : CONTAINS
164 : !!***
165 :
166 : !----------------------------------------------------------------------
167 :
168 : !!****f* m_eig2d/eigr2d_init
169 : !! NAME
170 : !! eigr2d_init
171 : !!
172 : !! FUNCTION
173 : !! This subroutine initializes the eigr2d_t structured datatype
174 : !!
175 : !! INPUTS
176 : !! mbands=maximum number of bands
177 : !! nkpt=number of k points
178 : !! nsppol=1 for unpolarized, 2 for spin-polarized
179 : !! natom=number of atoms
180 : !! eig2nkq(2,mband*nsppol,nkpt,3,natom,3,natom)=second-order derivative of the
181 : !! eigenenergies wrt phononic displacements
182 : !!
183 : !! OUTPUT
184 : !! eigr2d<eigr2d_t>=the eigr2d_t datatype
185 : !!
186 : !! SOURCE
187 :
188 0 : subroutine eigr2d_init(eig2nkq,eigr2d,mband,nsppol,nkpt,natom)
189 :
190 : !Arguments ------------------------------------
191 : !scalars
192 : integer,intent(in) ::mband,nsppol,nkpt,natom
193 : type(eigr2d_t),intent(out) :: eigr2d
194 : !arrays
195 : real(dp), intent(in) :: eig2nkq(2,mband*nsppol,nkpt,3,natom,3,natom)
196 :
197 : ! *************************************************************************
198 :
199 0 : eigr2d%mband = mband
200 0 : eigr2d%nsppol = nsppol
201 0 : eigr2d%nkpt = nkpt
202 0 : eigr2d%natom = natom
203 :
204 0 : ABI_MALLOC(eigr2d%eigr2d ,(2,mband*nsppol,nkpt,3,natom,3,natom))
205 0 : eigr2d%eigr2d=eig2nkq
206 :
207 0 : end subroutine eigr2d_init
208 : !!***
209 :
210 : !----------------------------------------------------------------------
211 :
212 : !!****f* m_eig2d/eigr2d_ncwrite
213 : !! NAME
214 : !! eigr2d_ncwrite
215 : !!
216 : !! FUNCTION
217 : !! Writes the content of a eigr2d_t object to a NETCDF file
218 : !! according to the ETSF-IO specifications.
219 : !!
220 : !! INPUTS
221 : !! ncid =NC file handle
222 : !!
223 : !! OUTPUT
224 : !!
225 : !! SOURCE
226 :
227 0 : subroutine eigr2d_ncwrite(eigr2d,iqpt,wtq,ncid)
228 :
229 : !Arguments ------------------------------------
230 : !scalars
231 : integer,intent(in) ::ncid
232 : real(dp),intent(in) :: iqpt(3),wtq
233 : type(eigr2d_t),intent(in) :: eigr2d
234 :
235 : !Local variables-------------------------------
236 : integer :: ncerr
237 : integer :: cplex,cart_dir,one_dim
238 : character(len=200) :: temp
239 : ! *************************************************************************
240 :
241 : ! ==============================================
242 : ! === Write the dimensions specified by ETSF ===
243 : ! ==============================================
244 0 : one_dim=1; cplex=2; cart_dir=3
245 :
246 : ncerr = nctk_def_dims(ncid, [&
247 : nctkdim_t('max_number_of_states', eigr2d%mband),&
248 : nctkdim_t('number_of_spins', eigr2d%nsppol),&
249 : nctkdim_t('number_of_kpoints', eigr2d%nkpt),&
250 : nctkdim_t('number_of_atoms', eigr2d%natom),&
251 : nctkdim_t('number_of_cartesian_directions', cart_dir),&
252 : nctkdim_t('current_one_dim', one_dim),&
253 : nctkdim_t('cplex', cplex),&
254 0 : nctkdim_t('product_mband_nsppol', eigr2d%mband*eigr2d%nsppol)], defmode=.True.)
255 0 : NCF_CHECK(ncerr)
256 :
257 : temp='cplex,product_mband_nsppol,number_of_kpoints,number_of_cartesian_directions,number_of_atoms,' //&
258 0 : 'number_of_cartesian_directions , number_of_atoms'
259 : ncerr = nctk_def_arrays(ncid, [&
260 : nctkarr_t('current_q_point', "dp", 'number_of_cartesian_directions'), &
261 : nctkarr_t('current_q_point_weight', "dp", 'current_one_dim'), &
262 0 : nctkarr_t('second_derivative_eigenenergies', "dp", temp )])
263 : ! nctkarr_t('second_derivative_eigenenergies', "dp",&
264 : ! &'cplex, product_mband_nsppol, number_of_kpoints, number_of_cartesian_directions, number_of_atoms,&
265 : ! &number_of_cartesian_directions, number_of_atoms')])
266 0 : NCF_CHECK(ncerr)
267 :
268 : ! Write data
269 0 : NCF_CHECK(nctk_set_datamode(ncid))
270 0 : NCF_CHECK(nf90_put_var(ncid, vid('current_q_point'), iqpt))
271 0 : NCF_CHECK(nf90_put_var(ncid, vid('current_q_point_weight'), wtq))
272 0 : NCF_CHECK(nf90_put_var(ncid, vid('second_derivative_eigenenergies'), eigr2d%eigr2d))
273 :
274 : contains
275 0 : integer function vid(vname)
276 : character(len=*),intent(in) :: vname
277 0 : vid = nctk_idname(ncid, vname)
278 : end function vid
279 :
280 : end subroutine eigr2d_ncwrite
281 : !!***
282 :
283 : !----------------------------------------------------------------------
284 :
285 : !!****f* m_eig2d/eigr2d_free
286 : !! NAME
287 : !! eigr2d_free
288 : !!
289 : !! FUNCTION
290 : !! Deallocates the components of the eigr2d_t structured datatype
291 : !!
292 : !! INPUTS
293 : !! eigr2d<eigr2d_t>=The data type to be deallocated.
294 : !!
295 : !! OUTPUT
296 : !! Deallocate the dynamic arrays in the ebands_t type.
297 : !! (only deallocate)
298 : !!
299 : !! SOURCE
300 :
301 0 : subroutine eigr2d_free(eigr2d)
302 :
303 : !Arguments ------------------------------------
304 : !scalars
305 : type(eigr2d_t),intent(inout) :: eigr2d
306 : ! *************************************************************************
307 : DBG_ENTER("COLL")
308 :
309 : !Deallocate all components of bstruct
310 0 : ABI_SFREE(eigr2d%eigr2d)
311 :
312 : DBG_EXIT("COLL")
313 :
314 0 : end subroutine eigr2d_free
315 : !!***
316 :
317 : !!****f* m_eig2d/fan_init
318 : !! NAME
319 : !! fan_init
320 : !!
321 : !! FUNCTION
322 : !! This subroutine initializes the fan_t structured datatype
323 : !!
324 : !! INPUTS
325 : !! mbands=maximum number of bands
326 : !! nkpt=number of k points
327 : !! nsppol=1 for unpolarized, 2 for spin-polarized
328 : !! natom=number of atoms
329 : !! fan2d(2*mband*nsppol,nkpt,3,natom,3,natom,mband*nsppol)=second-order derivative of the
330 : !! eigenenergies wrt phononic displacements
331 : !!
332 : !! OUTPUT
333 : !! fan2d<fan_t>=the fan_t datatype
334 : !!
335 : !! SIDE EFFECTS
336 : !!
337 : !! SOURCE
338 :
339 25 : subroutine fan_init(fan,fan2d,mband,nsppol,nkpt,natom)
340 :
341 : !Arguments ------------------------------------
342 : !scalars
343 : integer,intent(in) ::mband,nsppol,nkpt,natom
344 : type(fan_t),intent(out) :: fan2d
345 : !arrays
346 : real(dp), intent(in) :: fan(2*mband*nsppol,nkpt,3,natom,3,natom,mband)
347 : ! *************************************************************************
348 :
349 25 : fan2d%mband = mband
350 25 : fan2d%nsppol = nsppol
351 25 : fan2d%nkpt = nkpt
352 25 : fan2d%natom = natom
353 :
354 250 : ABI_MALLOC(fan2d%fan2d,(2*mband*nsppol,nkpt,3,natom,3,natom,mband))
355 1604402 : fan2d%fan2d=fan
356 :
357 25 : end subroutine fan_init
358 : !!***
359 :
360 : !!****f* m_eig2d/gkk_init
361 : !! NAME
362 : !! gkk_init
363 : !!
364 : !! FUNCTION
365 : !! This subroutine initializes the gkk_t structured datatype
366 : !!
367 : !! INPUTS
368 : !! mbands=maximum number of bands
369 : !! nkpt=number of k points
370 : !! nsppol=1 for unpolarized, 2 for spin-polarized
371 : !! natom=number of atoms
372 : !! gkk2d(2*mband*nsppol,nkpt,3,natom,mband*nsppol)=second-order derivative of the
373 : !! eigenenergies wrt phononic displacements
374 : !!
375 : !! OUTPUT
376 : !! gkk2d<gkk_t>=the gkk_t datatype
377 : !!
378 : !! SIDE EFFECTS
379 : !!
380 : !! SOURCE
381 :
382 122 : subroutine gkk_init(gkk,gkk2d,mband,nsppol,nkpt,natom,ncart)
383 :
384 : !Arguments ------------------------------------
385 : !scalars
386 : integer,intent(in) ::mband,nsppol,nkpt,natom,ncart
387 : type(gkk_t),intent(out) :: gkk2d
388 : !arrays
389 : real(dp), intent(in) :: gkk(2*mband*nsppol,nkpt,ncart,natom,mband)
390 : ! *************************************************************************
391 :
392 122 : gkk2d%mband = mband
393 122 : gkk2d%nsppol = nsppol
394 122 : gkk2d%nkpt = nkpt
395 122 : gkk2d%natom = natom
396 122 : gkk2d%ncart = ncart
397 :
398 976 : ABI_MALLOC(gkk2d%gkk2d,(2*mband*nsppol,nkpt,ncart,natom,mband))
399 2740315 : gkk2d%gkk2d=gkk
400 :
401 122 : end subroutine gkk_init
402 : !!***
403 :
404 : !----------------------------------------------------------------------
405 :
406 : !!****f* m_eig2d/fan_ncwrite
407 : !! NAME
408 : !! fan_ncwrite
409 : !!
410 : !! FUNCTION
411 : !! Writes the content of a fan_t object to a NETCDF file
412 : !! according to the ETSF-IO specifications.
413 : !!
414 : !! INPUTS
415 : !! ncid =NC file handle
416 : !!
417 : !! OUTPUT
418 : !!
419 : !! SOURCE
420 :
421 25 : subroutine fan_ncwrite(fan2d,iqpt,wtq,ncid)
422 :
423 : !Arguments ------------------------------------
424 : !scalars
425 : integer,intent(in) ::ncid
426 : real(dp),intent(in) :: iqpt(3),wtq
427 : type(fan_t),intent(in) :: fan2d
428 :
429 : !Local variables-------------------------------
430 : integer :: ncerr
431 : integer :: cplex,cart_dir,one_dim
432 : character(len=200) :: temp
433 : ! *************************************************************************
434 :
435 : ! ==============================================
436 : ! === Write the dimensions specified by ETSF ===
437 : ! ==============================================
438 25 : one_dim=1; cplex=2; cart_dir=3
439 :
440 : ncerr = nctk_def_dims(ncid, [&
441 : nctkdim_t('max_number_of_states',fan2d%mband),&
442 : nctkdim_t('number_of_spins',fan2d%nsppol),&
443 : nctkdim_t('number_of_kpoints',fan2d%nkpt),&
444 : nctkdim_t('number_of_atoms',fan2d%natom),&
445 : nctkdim_t('3_number_of_atoms',3*fan2d%natom),& ! TODO: not sure that variables can start with digits
446 : nctkdim_t('number_of_cartesian_directions',cart_dir),&
447 : nctkdim_t('current_one_dim',one_dim),&
448 : nctkdim_t('cplex',cplex),&
449 : nctkdim_t('product_mband_nsppol',fan2d%mband*fan2d%nsppol),&
450 : nctkdim_t('product_mband_nsppol2',fan2d%mband*fan2d%nsppol*2) &
451 275 : ], defmode=.True.)
452 25 : NCF_CHECK(ncerr)
453 :
454 : temp= 'product_mband_nsppol2, number_of_kpoints, number_of_cartesian_directions,' //&
455 25 : 'number_of_atoms, number_of_cartesian_directions, number_of_atoms, max_number_of_states'
456 : ncerr = nctk_def_arrays(ncid, [&
457 : nctkarr_t('current_q_point', "dp", 'number_of_cartesian_directions'),&
458 : nctkarr_t('current_q_point_weight', "dp", 'current_one_dim'),&
459 100 : nctkarr_t('second_derivative_eigenenergies_actif', "dp", temp )])
460 : ! nctkarr_t('second_derivative_eigenenergies_actif', "dp",&
461 : ! &'product_mband_nsppol2, number_of_kpoints, number_of_cartesian_directions,&
462 : ! &number_of_atoms, number_of_cartesian_directions, number_of_atoms, max_number_of_states')])
463 25 : NCF_CHECK(ncerr)
464 :
465 : ! Write data
466 25 : NCF_CHECK(nctk_set_datamode(ncid))
467 25 : NCF_CHECK(nf90_put_var(ncid, vid('current_q_point'), iqpt))
468 25 : NCF_CHECK(nf90_put_var(ncid, vid('current_q_point_weight'), wtq))
469 25 : NCF_CHECK(nf90_put_var(ncid, vid('second_derivative_eigenenergies_actif'), fan2d%fan2d))
470 :
471 : contains
472 75 : integer function vid(vname)
473 : character(len=*),intent(in) :: vname
474 75 : vid = nctk_idname(ncid, vname)
475 : end function vid
476 :
477 : end subroutine fan_ncwrite
478 : !!***
479 :
480 : !----------------------------------------------------------------------
481 :
482 : !!****f* m_eig2d/gkk_ncwrite
483 : !! NAME
484 : !! gkk_ncwrite
485 : !!
486 : !! FUNCTION
487 : !! Writes the content of a gkk_t object to a NETCDF file
488 : !! according to the ETSF-IO specifications.
489 : !!
490 : !! INPUTS
491 : !! ncid =NC file handle
492 : !!
493 : !! OUTPUT
494 : !!
495 : !! SOURCE
496 :
497 122 : subroutine gkk_ncwrite(gkk2d,iqpt,wtq,ncid)
498 :
499 : !Arguments ------------------------------------
500 : !scalars
501 : integer,intent(in) ::ncid
502 : real(dp), intent(in) :: iqpt(3),wtq
503 : type(gkk_t),intent(in) :: gkk2d
504 :
505 : !Local variables-------------------------------
506 : integer :: cplex,one_dim,ncerr,vid_
507 : ! *************************************************************************
508 :
509 : ! ==============================================
510 : ! === Write the dimensions specified by ETSF ===
511 : ! ==============================================
512 122 : one_dim=1; cplex=2
513 :
514 : ncerr = nctk_def_dims(ncid, [ &
515 : & nctkdim_t('max_number_of_states', gkk2d%mband), &
516 : & nctkdim_t('number_of_spins', gkk2d%nsppol), &
517 : & nctkdim_t('number_of_kpoints', gkk2d%nkpt), &
518 : & nctkdim_t('number_of_atoms_for_gkk', gkk2d%natom), &
519 : & nctkdim_t('3_number_of_atoms', 3*gkk2d%natom), &
520 : & nctkdim_t('number_of_cartesian_directions_for_gkk', gkk2d%ncart), &
521 : & nctkdim_t('current_one_dim', one_dim), &
522 : & nctkdim_t('cplex', cplex), &
523 : & nctkdim_t('product_mband_nsppol', gkk2d%mband*gkk2d%nsppol), &
524 : & nctkdim_t('product_mband_nsppol2', gkk2d%mband*gkk2d%nsppol*2) &
525 1342 : & ], defmode=.True.)
526 122 : NCF_CHECK(ncerr)
527 :
528 : !arrays
529 : ncerr = nctk_def_arrays(ncid, [&
530 : & nctkarr_t('current_q_point', "dp", "number_of_cartesian_directions"), &
531 : & nctkarr_t('current_q_point_weight', "dp", 'current_one_dim'), &
532 : & nctkarr_t('second_derivative_eigenenergies_actif', "dp", &
533 : & 'product_mband_nsppol2, number_of_kpoints, number_of_cartesian_directions_for_gkk,'// &
534 : & 'number_of_atoms_for_gkk, max_number_of_states') &
535 488 : & ])
536 122 : NCF_CHECK(ncerr)
537 :
538 122 : NCF_CHECK(nctk_set_datamode(ncid))
539 122 : vid_=vid('current_q_point')
540 122 : NCF_CHECK(nf90_put_var(ncid, vid_, iqpt))
541 122 : vid_=vid('current_q_point_weight')
542 122 : NCF_CHECK(nf90_put_var(ncid, vid_, wtq))
543 122 : vid_=vid('second_derivative_eigenenergies_actif')
544 122 : NCF_CHECK(nf90_put_var(ncid, vid_, gkk2d%gkk2d))
545 :
546 : contains
547 366 : integer function vid(vname)
548 : character(len=*),intent(in) :: vname
549 366 : vid = nctk_idname(ncid, vname)
550 : end function vid
551 :
552 : end subroutine gkk_ncwrite
553 : !!***
554 :
555 : !----------------------------------------------------------------------
556 :
557 : !!****f* m_eig2d/fan_free
558 : !! NAME
559 : !! fan_free
560 : !!
561 : !! FUNCTION
562 : !! Deallocates the components of the fan_t structured datatype
563 : !!
564 : !! INPUTS
565 : !! fan2d<fan_t>=The data type to be deallocated.
566 : !!
567 : !! OUTPUT
568 : !! Deallocate the dynamic arrays in the fan_t type.
569 : !! (only deallocate)
570 : !!
571 : !! SOURCE
572 :
573 54 : subroutine fan_free(fan2d)
574 :
575 : !Arguments ------------------------------------
576 : !scalars
577 : type(fan_t),intent(inout) :: fan2d
578 : ! *************************************************************************
579 : DBG_ENTER("COLL")
580 :
581 : !Deallocate all components of bstruct
582 :
583 54 : ABI_SFREE(fan2d%fan2d)
584 :
585 : DBG_EXIT("COLL")
586 :
587 0 : end subroutine fan_free
588 : !!***
589 :
590 : !----------------------------------------------------------------------
591 :
592 : !!****f* m_eig2d/gkk_free
593 : !! NAME
594 : !! gkk_free
595 : !!
596 : !! FUNCTION
597 : !! Deallocates the components of the gkk_t structured datatype
598 : !!
599 : !! INPUTS
600 : !! gkk2d<gkk_t>=The data type to be deallocated.
601 : !!
602 : !! OUTPUT
603 : !! Deallocate the dynamic arrays in the gkk_t type.
604 : !! (only deallocate)
605 : !!
606 : !! SOURCE
607 :
608 147 : subroutine gkk_free(gkk2d)
609 :
610 : !Arguments ------------------------------------
611 : !scalars
612 : type(gkk_t),intent(inout) :: gkk2d
613 : ! *************************************************************************
614 : DBG_ENTER("COLL")
615 :
616 : !Deallocate all components of bstruct
617 :
618 122 : ABI_SFREE(gkk2d%gkk2d)
619 :
620 : DBG_EXIT("COLL")
621 :
622 93 : end subroutine gkk_free
623 : !!***
624 :
625 : !!****f* ABINIT/eig2stern
626 : !! NAME
627 : !! eig2stern
628 : !!
629 : !! FUNCTION
630 : !! This routine calculates the second-order eigenvalues.
631 : !! The output eig2nkq is this quantity for the input k points.
632 : !!
633 : !! INPUTS
634 : !! bdeigrf = number of bands for which to calculate the second-order eigenvalues.
635 : !! clflg(3,mpert)= array on calculated perturbations for eig2rf.
636 : !! dim_eig2nkq = 1 if eig2nkq is to be computed.
637 : !! cg1_pert(2,mpw1*nspinor*mband*mk1mem*nsppol,3,mpert) = first-order wf in G
638 : !! space for each perturbation. The wavefunction is orthogonal to the
639 : !! active space.
640 : !! gh0c1_pert(2,mpw1*nspinor*mband*mk1mem*nsppol,3,mpert) = matrix containing the
641 : !! vector: <G|H(0)|psi(1)>, for each perturbation.
642 : !! gh1c_pert(2,mpw1*nspinor*mband*mk1mem*nsppol,3,mpert)) = matrix containing the
643 : !! vector: <G|H(1)|n,k>, for each perturbation. The wavefunction is
644 : !! orthogonal to the active space.
645 : !! eigbrd(2,mband*nsppol,nkpt,3,natom,3,natom) = broadening factors for the
646 : !! electronic eigenvalues (optional).
647 : !! eigen0(nkpt_rbz*mband*nsppol) = 0-order eigenvalues at all K-points:
648 : !! <k,n'|H(0)|k,n'> (hartree).
649 : !! eigenq(nkpt_rbz*mband*nsppol) = 0-order eigenvalues at all shifted K-points:
650 : !! <k+Q,n'|H(0)|k+Q,n'> (hartree).
651 : !! eigen1(nkpt_rbz*2*nsppol*mband**2,3,mpert) = matrix of first-order:
652 : !! <k+Q,n'|H(1)|k,n> (hartree) (calculated in dfpt_cgwf).
653 : !! eig2nkq(2,mband*nsppol,nkpt,3,natom,3,natom*dim_eig2nkq) = second derivatives of
654 : !! the electronic eigenvalues.
655 : !! elph2_imagden = imaginary part of the denominator of the sum-over-state expression
656 : !! for the electronic eigenenergy shift due to second-order electron-phonon
657 : !! interation.
658 : !! ieig2rf = integer for calculation type.
659 : !! indsym(4,nsym,natom) = indirect indexing array for atom labels
660 : !! (not used yet, but will be used with symmetries).
661 : !! istwfk_pert(nkpt_rbz,3,mpert) = integer for choice of storage of wavefunction at
662 : !! each k point for each perturbation.
663 : !! mband = maximum number of bands.
664 : !! mk1mem = maximum number of k points which can fit in memory (RF data);
665 : !! 0 if use disk.
666 : !! mpert = maximum number of perturbations.
667 : !! natom = number of atoms in the unit cell.
668 : !! npert = number of phonon perturbations, without taking into account directions:
669 : !! natom.
670 : !! nsym = number of symmetries (not used yet).
671 : !! mpi_enreg = information about MPI parallelization.
672 : !! mpw1 = maximum number of planewaves used to represent first-order wavefunctions.
673 : !! nkpt_rbz = number of k-points for each perturbation.
674 : !! npwar1(nkpt_rbz,mpert) = number of planewaves at k-point for first-order.
675 : !! nspinor = number of spinorial components of the wavefunctions.
676 : !! nsppol = 1 for unpolarized, 2 for spin-polarized.
677 : !! occ(mband*nkpt*nsppol)=occup number for each band (often 2) at each k point
678 : !! smdelta = integer controling the calculation of electron lifetimes.
679 : !! symq(4,2,nsym) = 1 if symmetry preserves present qpoint. From littlegroup_q (not used yet).
680 : !! symrec(3,3,nsym) = 3x3 matrices of the group symmetries (reciprocal space)
681 : !! (not used yet).
682 : !! symrel(3,3,nsym) = array containing the symmetries in real space (not used yet).
683 : !! timrev = 1 if time-reversal preserves the q wavevector; 0 otherwise
684 : !! (not in use yet).
685 : !! dtset = OPTIONAL, dataset structure containing the input variable of the
686 : !! calculation. This is required to use the k-interpolation routine.
687 : !! eigenq_fine(mband_fine,mkpt_fine,nsppol_fine) = OPTIONAL, 0-order eigenvalues
688 : !! at all shifted K-points: <k+Q,n'|H(0)|k+Q,n'> (hartree) of the
689 : !! fine grid. This information is read from the WF dense k-grid file.
690 : !! hdr_fine = OPTIONAL, header of the WF file of the fine k-point grid. This
691 : !! variable is required for the k-interpolation routine.
692 : !! hdr0 = OPTIONAL, header of the GS WF file of the corse k-point grid. This
693 : !! variable is required for the k-interpolation routine.
694 : !!
695 : !! OUTPUT
696 : !! eig2nkq(2,mband*nsppol,nkpt_rbz,3,npert,3,npert)= diagonal part of the
697 : !! second-order eigenvalues: E^{(2),diag}_{k,q,j}.
698 : !! eigbrd(2,mband*nsppol,nkpt_rbz,3,npert,3,npert)= OPTIONAL, array containing the
699 : !! electron lifetimes.
700 : !!
701 : !! SOURCE
702 :
703 168 : subroutine eig2stern(dtfil,occ,bdeigrf,clflg,cg1_pert,dim_eig2nkq,dim_eig2rf,eigen0,eigenq,&
704 168 : & eigen1,eig2nkq,elph2_imagden,esmear,gh0c1_pert,gh1c_pert,ieig2rf,istwfk_pert,&
705 84 : & mband,mk1mem,mpert,npert,mpi_enreg,mpw1,nkpt_rbz,npwar1,nspinor,nsppol,smdelta,&
706 82 : & dtset,xred,pawtab,psps,eigbrd,eigenq_fine,hdr_fine,hdr0)
707 :
708 : !Arguments ------------------------------------
709 : type(datafiles_type), intent(in) :: dtfil
710 : type(MPI_type), intent(inout) :: mpi_enreg
711 : !scalars
712 : integer,intent(in) :: bdeigrf,dim_eig2nkq,dim_eig2rf,ieig2rf,mband,mk1mem,mpert,mpw1,nkpt_rbz
713 : integer,intent(in) :: npert,nspinor,nsppol,smdelta
714 : integer,intent(in) :: clflg(3,mpert)
715 : integer,intent(in) :: istwfk_pert(nkpt_rbz,3,mpert)
716 : integer,intent(in) :: npwar1(nkpt_rbz,mpert)
717 : !arrays
718 : type(dataset_type), intent(in) :: dtset
719 : real(dp),intent(in) :: elph2_imagden,esmear
720 : real(dp),intent(in) :: cg1_pert(2,mpw1*nspinor*mband*mk1mem*nsppol*dim_eig2rf,3,mpert)
721 : real(dp),intent(in) :: gh0c1_pert(2,mpw1*nspinor*mband*mk1mem*nsppol*dim_eig2rf,3,mpert)
722 : real(dp),intent(in) :: gh1c_pert(2,mpw1*nspinor*mband*mk1mem*nsppol*dim_eig2rf,3,mpert)
723 : real(dp),intent(inout) :: eigen0(nkpt_rbz*mband*nsppol)
724 : real(dp),intent(inout) :: xred(3,dtset%natom)
725 : real(dp),intent(in) :: eigen1(nkpt_rbz*2*nsppol*mband**2,3,mpert)
726 : real(dp),intent(inout) :: eigenq(nkpt_rbz*mband*nsppol)
727 : real(dp),intent(out) :: eig2nkq(2,mband*nsppol,nkpt_rbz,3,npert,3,npert*dim_eig2nkq)
728 : real(dp),intent(out),optional :: eigbrd(2,mband*nsppol,nkpt_rbz,3,npert,3,npert)
729 : real(dp),intent(in),pointer,optional :: eigenq_fine(:,:,:)
730 : real(dp), intent(in) :: occ(mband*nkpt_rbz*nsppol)
731 : type(pseudopotential_type), intent(inout) :: psps
732 : type(pawtab_type),intent(inout) :: pawtab(psps%ntypat*psps%usepaw)
733 : type(hdr_type),intent(in),optional :: hdr_fine,hdr0
734 :
735 : !Local variables-------------------------------
736 : !tolerance for non degenerated levels
737 : !scalars
738 : integer :: band2tot_index,band_index,bandtot_index,iband,icg2,idir1,idir2
739 : integer :: ikpt,ipert1,ipert2,isppol,istwf_k,jband,npw1_k,nkpt_sub,ikpt2
740 : !integer :: ipw
741 : integer :: master,me,spaceworld,ierr,mband_mem,mpert_
742 : !real(dp),parameter :: etol=1.0d-3
743 : real(dp),parameter :: etol=1.0d-6
744 : !real(dp),parameter :: etol=zero
745 : real(dp) :: ar,ai,deltae,den,dot2i,dot2r,dot3i,dot3r,doti,dotr,eig1_i1,eig1_i2
746 : real(dp) :: eig1_r1,eig1_r2,eig2_diai,den_av
747 : real(dp) :: wgt_int
748 : real(dp) :: eig2_diar,eigbrd_i,eigbrd_r
749 : character(len=500) :: message
750 : character(len=500) :: msg
751 : character(len=fnlen) :: dscrpt
752 : !DBSP
753 : ! character(len=300000) :: message2
754 : !END
755 : logical :: test_do_band
756 : !arrays
757 84 : integer, allocatable :: nband_rbz(:),icg2_rbz(:,:)
758 84 : integer, allocatable :: blkflg_save(:,:,:,:)
759 84 : integer, allocatable :: flg(:,:,:,:)
760 84 : real(dp),allocatable :: eig2nkq_tmp(:,:,:,:,:,:,:)
761 84 : integer, pointer :: kpt_fine_sub(:)
762 : real(dp) :: tsec(2)
763 168 : real(dp),allocatable :: cwavef(:,:),cwavef2(:,:),center(:),eigen0tmp(:),eigenqtmp(:)
764 252 : real(dp) :: eigen(mband*nsppol),eigen_prime(mband*nsppol)
765 84 : real(dp),allocatable :: gh(:,:),gh1(:,:),ghc(:,:)
766 84 : real(dp),allocatable :: smdfun(:,:)
767 84 : real(dp),pointer :: wgt_sub(:)
768 84 : type(ddb_type) :: ddb
769 4368 : type(ddb_hdr_type) :: ddb_hdr
770 :
771 : ! *********************************************************************
772 :
773 : !Init parallelism
774 84 : master =0
775 84 : spaceworld=mpi_enreg%comm_cell
776 84 : me=mpi_enreg%me_kpt
777 :
778 : !Init interpolation method
779 84 : if(present(eigenq_fine))then
780 3 : ABI_MALLOC(center,(3))
781 : end if
782 :
783 84 : call timab(148,1,tsec)
784 :
785 84 : if(nsppol==2)then
786 1 : message = 'nsppol=2 is still under development. Be careful when using it ...'
787 1 : ABI_COMMENT(message)
788 : end if
789 :
790 84 : band2tot_index =0
791 84 : bandtot_index=0
792 84 : band_index=0
793 :
794 : !Add scissor shift to eigenenergies
795 84 : if (dtset%dfpt_sciss > tol6 ) then
796 : write(msg,'(a,f7.3,2a)')&
797 0 : & ' A scissor operator of ',dtset%dfpt_sciss*Ha_eV,' [eV] has been applied to the eigenenergies',ch10
798 0 : call wrtout(std_out,msg,'COLL')
799 0 : call wrtout(ab_out,msg,'COLL')
800 0 : ABI_MALLOC(eigen0tmp,(nkpt_rbz*mband*nsppol))
801 0 : ABI_MALLOC(eigenqtmp,(nkpt_rbz*mband*nsppol))
802 0 : eigen0tmp = eigen0(:)
803 0 : eigenqtmp = eigenq(:)
804 0 : eigen0 = zero
805 0 : eigenq = zero
806 : end if
807 :
808 84 : if(ieig2rf > 0) then
809 1723672 : eig2nkq(:,:,:,:,:,:,:) = zero
810 : end if
811 84 : if(present(eigbrd))then
812 1686694 : eigbrd(:,:,:,:,:,:,:) = zero
813 : end if
814 :
815 : if(xmpi_paral==1) then
816 420 : ABI_MALLOC(mpi_enreg%proc_distrb,(nkpt_rbz,mband,nsppol))
817 252 : ABI_MALLOC(nband_rbz,(nkpt_rbz*nsppol))
818 84 : ABI_SFREE(mpi_enreg%my_kpttab)
819 252 : ABI_MALLOC(mpi_enreg%my_kpttab,(nkpt_rbz))
820 : ! Assume the number of bands is the same for all k points.
821 1569 : nband_rbz(:)=mband
822 84 : call distrb2(mband,mband_mem,nband_rbz,nkpt_rbz,mpi_enreg%nproc_cell,nsppol,mpi_enreg)
823 : end if
824 :
825 84 : icg2=0
826 84 : ipert1=1 ! Suppose that the situation is the same for all perturbations
827 336 : ABI_MALLOC(icg2_rbz,(nkpt_rbz,nsppol))
828 169 : do isppol=1,nsppol
829 1654 : do ikpt=1,nkpt_rbz
830 1485 : icg2_rbz(ikpt,isppol)=icg2
831 1485 : if(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,mband,isppol,me)) cycle
832 1570 : icg2 = icg2 + npwar1(ikpt,ipert1)*nspinor*mband
833 : end do
834 : end do
835 :
836 169 : do isppol=1,nsppol
837 1570 : do ikpt =1,nkpt_rbz
838 :
839 1485 : if(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,mband,isppol,me)) then
840 224 : band2tot_index = band2tot_index + 2*mband**2
841 224 : bandtot_index = bandtot_index + mband
842 224 : cycle
843 : end if
844 :
845 1261 : if(present(eigenq_fine))then
846 24 : write(std_out,*) 'Start of the energy denominator interpolation method.'
847 24 : nkpt_sub = 0
848 : ! center is the k+q point around which we will average the kpt_fine
849 120 : center = hdr0%kptns(:,ikpt)+ dtset%qptn(:)
850 :
851 24 : call kptfine_av(center,dtset%qptrlatt,hdr_fine%kptns,hdr_fine%nkpt,kpt_fine_sub,nkpt_sub,wgt_sub)
852 : write(std_out,'(a,3f8.4,a,i3)') 'Number of k-points of the fine grid &
853 24 : & around the k+Q point ',center,' is:',nkpt_sub
854 672 : write(std_out,'(a,f10.5)') 'The sum of the weights of the k-points is: ',SUM(wgt_sub)
855 : end if
856 :
857 : ! Add scissor shift to eigenenergies
858 1261 : if (dtset%dfpt_sciss > tol6 ) then
859 0 : do iband=1,mband
860 0 : if (occ(iband+bandtot_index) < tol6) then
861 0 : eigen0(iband+bandtot_index) = eigen0tmp(iband+bandtot_index) + dtset%dfpt_sciss
862 0 : eigenq(iband+bandtot_index) = eigenqtmp(iband+bandtot_index) + dtset%dfpt_sciss
863 : else
864 0 : eigen0(iband+bandtot_index) = eigen0tmp(iband+bandtot_index)
865 0 : eigenq(iband+bandtot_index) = eigenqtmp(iband+bandtot_index)
866 : end if
867 : end do
868 : end if
869 :
870 :
871 1261 : if(smdelta >0) then !broadening
872 1197 : if(.not.allocated(smdfun)) then
873 1443 : ABI_MALLOC(smdfun,(mband,mband))
874 : end if
875 159073 : smdfun(:,:) = zero
876 14033 : do iband=1,mband
877 12836 : eigen(iband) = eigen0(iband+bandtot_index)
878 14033 : eigen_prime(iband) =eigenq(iband+bandtot_index)
879 : end do
880 1197 : if(esmear>tol6) then
881 1197 : call smeared_delta(eigen,eigen_prime,esmear,mband,smdelta,smdfun)
882 : end if
883 : end if
884 1261 : icg2=icg2_rbz(ikpt,isppol)
885 :
886 1261 : ipert1=1 ! Suppose all perturbations lead to the same number of planewaves
887 1261 : npw1_k = npwar1(ikpt,ipert1)
888 3783 : ABI_MALLOC(cwavef,(2,npw1_k*nspinor))
889 2522 : ABI_MALLOC(cwavef2,(2,npw1_k*nspinor))
890 2522 : ABI_MALLOC(gh,(2,npw1_k*nspinor))
891 2522 : ABI_MALLOC(gh1,(2,npw1_k*nspinor))
892 2522 : ABI_MALLOC(ghc,(2,npw1_k*nspinor))
893 :
894 14397 : do iband=1,bdeigrf
895 :
896 : ! If the k point and band belong to me, compute the contribution
897 13136 : test_do_band=.true.
898 13136 : if(mpi_enreg%proc_distrb(ikpt,iband,isppol)/=me)test_do_band=.false.
899 :
900 1261 : if(test_do_band)then
901 :
902 39344 : do ipert1=1,npert
903 :
904 117968 : do idir1=1,3
905 78624 : if(clflg(idir1,ipert1)==0)cycle
906 78584 : istwf_k = istwfk_pert(ikpt,idir1,ipert1)
907 :
908 265800 : do ipert2=1,npert
909 722656 : do idir2=1,3
910 483024 : if(clflg(idir2,ipert2)==0)cycle
911 :
912 : eig2_diar = zero ; eig2_diai = zero ; eigbrd_r = zero ; eigbrd_i = zero
913 :
914 6256096 : do jband=1,mband
915 5773152 : eig1_r1 = eigen1(2*jband-1+(iband-1)*2*mband+band2tot_index,idir1,ipert1)
916 5773152 : eig1_r2 = eigen1(2*jband-1+(iband-1)*2*mband+band2tot_index,idir2,ipert2)
917 5773152 : eig1_i1 = eigen1(2*jband+(iband-1)*2*mband+band2tot_index,idir1,ipert1)
918 5773152 : eig1_i2 = - eigen1(2*jband+(iband-1)*2*mband+band2tot_index,idir2,ipert2) !the negative sign is from the CC
919 : ! If no interpolation, fallback on to the previous
920 : ! implementation
921 5773152 : if(.not. present(eigenq_fine))then
922 5717856 : deltae=eigenq(jband+bandtot_index)-eigen0(iband+bandtot_index)
923 : end if
924 5773152 : ar=eig1_r1*eig1_r2-eig1_i1*eig1_i2
925 5773152 : ai=eig1_r1*eig1_i2+eig1_i1*eig1_r2
926 :
927 : ! Sum over all active space to retrieve the diagonal gauge
928 5773152 : if(ieig2rf == 1 .or. ieig2rf ==2 ) then
929 : ! if(abs(deltae)>etol) then ! This is commented because
930 : ! there is no problem with divergencies with elph2_imag != 0
931 1036992 : if( present(eigenq_fine))then
932 55296 : den_av = zero
933 55296 : wgt_int = zero
934 1548288 : do ikpt2=1,nkpt_sub
935 : deltae=eigenq_fine(jband,kpt_fine_sub(ikpt2),1)&
936 1492992 : & -eigen0(iband+bandtot_index)
937 1492992 : den_av = den_av-(wgt_sub(ikpt2)*deltae)/(deltae**2+elph2_imagden**2)
938 1548288 : wgt_int = wgt_int+wgt_sub(ikpt2)
939 : end do
940 55296 : den = den_av/wgt_int
941 : else
942 981696 : if(abs(elph2_imagden) < etol) then
943 723648 : if(abs(deltae)>etol) then
944 548132 : den=-one/(deltae**2+elph2_imagden**2)
945 : else
946 : den= zero
947 : end if
948 : else
949 258048 : den=-one/(deltae**2+elph2_imagden**2)
950 : end if
951 : end if
952 :
953 : ! The following should be the most general implementation of the presence of elph2_imagden
954 : ! eig2_diar=eig2_diar+(ar*deltae+ai*elph2_imagden)*den
955 : ! eig2_diai=eig2_diai+(ai*deltae-ar*elph2_imagden)*den
956 : ! This gives back the implementation without elph2_imagden
957 : ! eig2_diar=eig2_diar+ar*deltae*den
958 : ! eig2_diai=eig2_diai+ai*deltae*den
959 : ! This is what Samuel had implemented
960 : ! eig2_diar=eig2_diar+ar*deltae*den
961 : ! eig2_diai=eig2_diai+ai*elph2_imagden*den
962 : ! Other possibility : throw away the broadening part, that is actually treated separately.
963 1036992 : if( present(eigenq_fine))then
964 55296 : eig2_diar=eig2_diar+ar*den
965 55296 : eig2_diai=eig2_diai+ai*den
966 : else
967 981696 : eig2_diar=eig2_diar+ar*deltae*den
968 981696 : eig2_diai=eig2_diai+ai*deltae*den
969 : !DBSP
970 : ! if (iband+band_index==2 .and. ikpt==1 .and. idir1==1 .and. ipert1==1 .and. idir2==1 .and. ipert2==1) then
971 : ! write(message2,*) 'eig2_diar1=',eig2_diar,' ar=',ar,' deltae=',deltae,' den=',den
972 : ! call wrtout(std_out,message2,'PERS')
973 : ! endif
974 : !END
975 :
976 : end if
977 : end if ! ieig2rf==1 or 2
978 :
979 6256096 : if(present(eigbrd))then
980 5715552 : if(smdelta >0) then !broadening
981 5715552 : eigbrd_r = eigbrd_r + ar*smdfun(iband,jband)
982 5715552 : eigbrd_i = eigbrd_i + ai*smdfun(iband,jband)
983 : end if
984 : end if
985 :
986 : end do !jband
987 :
988 : ! Add the contribution of non-active bands, if DFPT calculation (= Sternheimer)
989 482944 : if(ieig2rf == 1 .or. ieig2rf ==3 .or. ieig2rf ==4 .or. ieig2rf==5 ) then
990 : ! if(ieig2rf == 1 ) then
991 :
992 482904 : dotr=zero ; doti=zero
993 482904 : dot2r=zero ; dot2i=zero
994 482904 : dot3r=zero ; dot3i=zero
995 :
996 :
997 923921136 : cwavef(:,:) = cg1_pert(:,1+(iband-1)*npw1_k*nspinor+icg2:iband*npw1_k*nspinor+icg2,idir2,ipert2)
998 923921136 : cwavef2(:,:)= cg1_pert(:,1+(iband-1)*npw1_k*nspinor+icg2:iband*npw1_k*nspinor+icg2,idir1,ipert1)
999 923921136 : gh1(:,:) = gh1c_pert(:,1+(iband-1)*npw1_k*nspinor+icg2:iband*npw1_k*nspinor+icg2,idir1,ipert1)
1000 923921136 : gh(:,:) = gh1c_pert(:,1+(iband-1)*npw1_k*nspinor+icg2:iband*npw1_k*nspinor+icg2,idir2,ipert2)
1001 923921136 : ghc(:,:) = gh0c1_pert(:,1+(iband-1)*npw1_k*nspinor+icg2:iband*npw1_k*nspinor+icg2,idir1,ipert1)
1002 :
1003 : ! The first two dotprod corresponds to: <Psi(1)nkq|H(1)k+q,k|Psi(0)nk> and <Psi(0)nk|H(1)k,k+q|Psi(1)nkq>
1004 : ! They are calculated using wavefunctions <Psi(1)| that are orthogonal to the active space.
1005 482904 : call dotprod_g(dotr,doti,istwf_k,npw1_k*nspinor,2,cwavef,gh1,mpi_enreg%me_g0,mpi_enreg%comm_spinorfft)
1006 482904 : call dotprod_g(dot2r,dot2i,istwf_k,npw1_k*nspinor,2,gh,cwavef2,mpi_enreg%me_g0,mpi_enreg%comm_spinorfft)
1007 :
1008 : ! This dotprod corresponds to : <Psi(1)nkq|H(0)k+q- E(0)nk|Psi(1)nkq>
1009 : ! It is calculated using wavefunctions that are orthogonal to the active space.
1010 : ! Should work for metals. (But adiabatic approximation is bad in this case...)
1011 482904 : call dotprod_g(dot3r,dot3i,istwf_k,npw1_k*nspinor,2,cwavef,ghc,mpi_enreg%me_g0,mpi_enreg%comm_spinorfft)
1012 :
1013 482904 : eig2_diar= eig2_diar + dotr + dot2r + dot3r
1014 482904 : eig2_diai= eig2_diai + doti + dot2i + dot3i
1015 :
1016 : end if
1017 :
1018 : ! Store the contribution
1019 482944 : if(ieig2rf > 0) then
1020 482944 : eig2nkq(1,iband+band_index,ikpt,idir1,ipert1,idir2,ipert2) = eig2_diar
1021 482944 : eig2nkq(2,iband+band_index,ikpt,idir1,ipert1,idir2,ipert2) = eig2_diai
1022 : end if
1023 :
1024 643952 : if(present(eigbrd))then
1025 471424 : if(smdelta >0) then !broadening
1026 471424 : eigbrd(1,iband+band_index,ikpt,idir1,ipert1,idir2,ipert2) = eigbrd_r
1027 471424 : eigbrd(2,iband+band_index,ikpt,idir1,ipert1,idir2,ipert2) = eigbrd_i
1028 : end if
1029 : end if
1030 :
1031 : end do !idir2
1032 : end do !ipert2
1033 : end do !idir1
1034 : end do !ipert1
1035 :
1036 : end if ! Selection of processor
1037 :
1038 : end do !iband
1039 :
1040 1261 : ABI_FREE(cwavef)
1041 1261 : ABI_FREE(cwavef2)
1042 1261 : ABI_FREE(gh)
1043 1261 : ABI_FREE(gh1)
1044 1261 : ABI_FREE(ghc)
1045 1261 : band2tot_index = band2tot_index + 2*mband**2
1046 1261 : bandtot_index = bandtot_index + mband
1047 :
1048 1346 : if(present(eigenq_fine))then
1049 24 : ABI_FREE(kpt_fine_sub) ! Deallocate the variable
1050 24 : ABI_FREE(wgt_sub)
1051 : end if
1052 :
1053 : end do !ikpt
1054 169 : band_index = band_index + mband
1055 : end do !isppol
1056 :
1057 : !Accumulate eig2nkq and/or eigbrd
1058 : if(xmpi_paral==1) then
1059 84 : if(ieig2rf == 1 .or. ieig2rf == 2) then
1060 30 : call xmpi_sum(eig2nkq,spaceworld,ierr)
1061 30 : if (dtset%dfpt_sciss > tol6 ) then
1062 0 : call xmpi_sum(eigen0,spaceworld,ierr)
1063 0 : call xmpi_sum(eigenq,spaceworld,ierr)
1064 : end if
1065 : end if
1066 84 : if(present(eigbrd) .and. (ieig2rf == 1 .or. ieig2rf == 2))then
1067 28 : if(smdelta >0) then
1068 28 : call xmpi_sum(eigbrd,spaceworld,ierr)
1069 : end if
1070 : end if
1071 84 : ABI_FREE(nband_rbz)
1072 84 : ABI_FREE(mpi_enreg%proc_distrb)
1073 84 : ABI_FREE(mpi_enreg%my_kpttab)
1074 : end if
1075 :
1076 84 : if(ieig2rf==1 .or. ieig2rf==2 ) then
1077 30 : write(ab_out,'(a)')' Components of second-order derivatives of the electronic energy, EIGR2D, in Ha unit.'
1078 30 : write(ab_out,'(a)')' For automatic tests, printing the matrix for the first k-point, first band, first atom.'
1079 120 : do idir1=1,3
1080 390 : do idir2=1,3
1081 270 : ar=eig2nkq(1,1,1,idir1,1,idir2,1) ; if(abs(ar)<tol10)ar=zero
1082 270 : ai=eig2nkq(2,1,1,idir1,1,idir2,1) ; if(abs(ai)<tol10)ai=zero
1083 360 : write (ab_out,'(4i4,2es20.10)') idir1,1,idir2,1,ar,ai
1084 : end do ! idir2
1085 : end do ! idir1
1086 : end if
1087 84 : if(ieig2rf > 2) then
1088 54 : write(ab_out,'(a)')' Components of the Sternheimer part only of the second-order derivatives of the electronic energy,'
1089 54 : write(ab_out,'(a)')' EIGR2D, in Ha unit.'
1090 54 : write(ab_out,'(a)')' For automatic tests, printing the matrix for the first k-point, first band, first atom.'
1091 216 : do idir1=1,3
1092 702 : do idir2=1,3
1093 486 : ar=eig2nkq(1,1,1,idir1,1,idir2,1) ; if(abs(ar)<tol10)ar=zero
1094 486 : ai=eig2nkq(2,1,1,idir1,1,idir2,1) ; if(abs(ai)<tol10)ai=zero
1095 648 : write (ab_out,'(4i4,2es20.10)') idir1,1,idir2,1,ar,ai
1096 : end do ! idir2
1097 : end do ! idir1
1098 : end if
1099 84 : if(present(eigbrd))then
1100 82 : if(smdelta >0) then !broadening
1101 82 : write(ab_out,'(a)')' '
1102 82 : write(ab_out,'(a)')' Components of second-order derivatives of the electronic energy, EIGI2D.'
1103 82 : write(ab_out,'(a)')' For automatic tests, printing the matrix for the first k-point, first band, first atom.'
1104 328 : do idir1=1,3
1105 1066 : do idir2=1,3
1106 738 : ar=eigbrd(1,1,1,idir1,1,idir2,1) ; if(abs(ar)<tol10)ar=zero
1107 738 : ai=eigbrd(2,1,1,idir1,1,idir2,1) ; if(abs(ai)<tol10)ai=zero
1108 984 : write (ab_out,'(4i4,2es20.10)') idir1,1,idir2,1,ar,ai
1109 : end do
1110 : end do !nband
1111 : end if
1112 : end if
1113 :
1114 : ! -------------------------
1115 : ! Output d2eig data to file
1116 : ! -------------------------
1117 : ! Write the second-order derivatives in reduced coordinate of the electronic energy, EIG2D to file, in Hartree unit.
1118 : ! When ieig2rf==1 or 2 --> EIG(R/I)2D contains the total second-order derivatives
1119 : ! When ieig2rf > 2 --> EIG(R/I)2D contains the Sternheimer part only.
1120 : !
1121 84 : if (dtset%ieig2rf==1 .or. dtset%ieig2rf==2) then
1122 : ! SP: For this case, the eig2nkq contains the active space and Sternheimer part
1123 : !
1124 : ! GA: Here, mpert needs to be replaced by natom
1125 : ! but why is mpert larger than natom in the first place?
1126 30 : mpert_ = dtset%natom
1127 :
1128 : ! Initialize perturbation flags
1129 : ! GA: At the moment, they are all set to one
1130 : ! Instead, they should be used to save individual perturbations
1131 : ! to separate files and merge them after the loop.
1132 150 : ABI_MALLOC(blkflg_save,(3,mpert_,3,mpert_))
1133 2134 : blkflg_save = one
1134 :
1135 : ! Initialize ddb object
1136 : call ddb%init(dtset, 1, mpert_, &
1137 : mband=bdeigrf,&
1138 : nkpt=nkpt_rbz,&
1139 : kpt=dtset%kptns(1:3,1:nkpt_rbz),&
1140 30 : with_d2eig=.true.)
1141 :
1142 : ! Create the ddb header
1143 30 : dscrpt=' Note : temporary (transfer) database '
1144 : call ddb_hdr%init(dtset,psps,pawtab,dscrpt,1,&
1145 : mpert=mpert_,&
1146 : xred=xred,occ=occ,&
1147 : mband=bdeigrf / dtset%nsppol,&
1148 : nkpt=nkpt_rbz,&
1149 30 : kpt=dtset%kptns(:,1:nkpt_rbz))
1150 :
1151 : ! Set d2eig data
1152 30 : call ddb%set_qpt(1, dtset%qptn)
1153 30 : call ddb%set_d2eig_reshape(1, eig2nkq, blkflg_save)
1154 :
1155 30 : call ddb_hdr%set_typ(ddb%nblok, ddb%typ)
1156 :
1157 : ! Open the file and write header
1158 30 : call ddb_hdr%open_write(dtfil%fnameabo_eigr2d, with_psps=1, comm=mpi_enreg%comm_world)
1159 :
1160 : ! Write d2eig data block
1161 30 : call ddb%write_d2eig(ddb_hdr, 1, comm=mpi_enreg%comm_world)
1162 :
1163 : ! close and free memory
1164 30 : call ddb_hdr%close()
1165 30 : call ddb_hdr%free()
1166 30 : call ddb%free()
1167 :
1168 30 : if (smdelta>0) then
1169 : ! write out _EIGI2D file
1170 :
1171 : call ddb%init(dtset, 1, mpert_, &
1172 : mband=bdeigrf,&
1173 : nkpt=nkpt_rbz,&
1174 : kpt=dtset%kptns(:,1:nkpt_rbz),&
1175 28 : with_d2eig=.true.)
1176 :
1177 : ! Create the ddb header
1178 28 : dscrpt=' Note : temporary (transfer) database '
1179 : call ddb_hdr%init(dtset,psps,pawtab,dscrpt,1,&
1180 : mpert=mpert_,&
1181 : xred=xred,occ=occ,&
1182 : mband=bdeigrf / dtset%nsppol,&
1183 : nkpt=nkpt_rbz,&
1184 28 : kpt=dtset%kptns(:,1:nkpt_rbz))
1185 :
1186 28 : call ddb%set_qpt(1, dtset%qptn)
1187 28 : call ddb%set_d2eig_reshape(1, eigbrd, blkflg_save, blktyp=BLKTYP_d2eig_im)
1188 :
1189 28 : call ddb_hdr%set_typ(ddb%nblok, ddb%typ)
1190 :
1191 28 : call ddb_hdr%open_write(dtfil%fnameabo_eigi2d, with_psps=1,comm=mpi_enreg%comm_world)
1192 28 : call ddb%write_d2eig(ddb_hdr, 1, comm=mpi_enreg%comm_world)
1193 :
1194 28 : call ddb_hdr%close()
1195 28 : call ddb_hdr%free()
1196 28 : call ddb%free()
1197 :
1198 : end if !smdelta
1199 :
1200 30 : ABI_FREE(blkflg_save)
1201 :
1202 : end if !ieig2rf==1.or.ieig2rf==2
1203 : !
1204 84 : if (ieig2rf == 3 .or. ieig2rf == 4 .or. ieig2rf == 5) then
1205 : ! SP: these cases are used for non-adiabatic quantities. In this routine we do not have access
1206 : ! to the phonon frequency. Therefore eig2nkq only contains the Sternheimer part.
1207 : !
1208 :
1209 54 : mpert_ = dtset%natom
1210 :
1211 : ! Initialize perturbation flags
1212 270 : ABI_MALLOC(flg,(3,mpert_,3,mpert_))
1213 3078 : flg = one
1214 :
1215 : ! Initialize ddb object
1216 : call ddb%init(dtset, 1, mpert_, &
1217 : mband=bdeigrf,&
1218 : nkpt=nkpt_rbz,&
1219 : kpt=dtset%kptns(1:3,1:nkpt_rbz),&
1220 54 : with_d2eig=.true.)
1221 :
1222 : ! Create the ddb header
1223 54 : dscrpt=' Note : temporary (transfer) database '
1224 : call ddb_hdr%init(dtset,psps,pawtab,dscrpt,1,&
1225 : mpert=mpert_,&
1226 : xred=xred,occ=occ,&
1227 : mband=bdeigrf / dtset%nsppol,&
1228 : nkpt=nkpt_rbz,&
1229 54 : kpt=dtset%kptns(:,1:nkpt_rbz))
1230 :
1231 : ! Set d2eig data
1232 54 : call ddb%set_qpt(1, dtset%qptn)
1233 :
1234 : ! Copy eig2nkq because it needs to be summed before writing the file,
1235 : ! yet there are further processing on the eig2nkq array occuring
1236 : ! after this function
1237 432 : ABI_MALLOC(eig2nkq_tmp,(2,mband*nsppol,nkpt_rbz,3,npert,3,npert*dim_eig2nkq))
1238 1347408 : eig2nkq_tmp = eig2nkq
1239 : if(xmpi_paral==1) then
1240 54 : call xmpi_sum(eig2nkq_tmp,spaceworld,ierr)
1241 : end if
1242 :
1243 54 : call ddb%set_d2eig_reshape(1, eig2nkq_tmp, flg)
1244 :
1245 54 : ABI_FREE(eig2nkq_tmp)
1246 :
1247 : ! Open the file and write header
1248 54 : call ddb_hdr%set_typ(ddb%nblok, ddb%typ)
1249 :
1250 : ! Write d2eig data block
1251 54 : call ddb_hdr%open_write(dtfil%fnameabo_eigr2d, with_psps=1, comm=mpi_enreg%comm_world)
1252 54 : call ddb%write_d2eig(ddb_hdr, 1, comm=mpi_enreg%comm_world)
1253 :
1254 : ! close and free memory
1255 54 : call ddb_hdr%close()
1256 54 : call ddb_hdr%free()
1257 54 : call ddb%free()
1258 : !
1259 108 : ABI_FREE(flg)
1260 : !
1261 : end if
1262 :
1263 84 : if(allocated(smdfun)) then
1264 82 : ABI_FREE(smdfun)
1265 : end if
1266 84 : ABI_FREE(icg2_rbz)
1267 84 : if(present(eigenq_fine))then
1268 3 : ABI_FREE(center)
1269 : end if
1270 84 : if (dtset%dfpt_sciss > tol6 ) then
1271 0 : ABI_FREE(eigen0tmp)
1272 0 : ABI_FREE(eigenqtmp)
1273 : end if
1274 :
1275 84 : call timab(148,2,tsec)
1276 :
1277 252 : end subroutine eig2stern
1278 : !!***
1279 :
1280 : !!****f* m_eig2d/eig2tot
1281 : !! NAME
1282 : !! eig2tot
1283 : !!
1284 : !! FUNCTION
1285 : !! This routine calculates the second-order eigenvalues.
1286 : !! The output eig2nkq is this quantity for the input k points.
1287 : !!
1288 : !! INPUTS
1289 : !! bdeigrf = number of bands for which to calculate the second-order eigenvalues.
1290 : !! clflg(3,mpert)= array on calculated perturbations for eig2rf.
1291 : !! dim_eig2nkq = 1 if eig2nkq is to be computed.
1292 : !! eigbrd(2,mband*nsppol,nkpt,3,natom,3,natom) = broadening factors for the
1293 : !! electronic eigenvalues (optional).
1294 : !! eigen0(nkpt_rbz*mband*nsppol) = 0-order eigenvalues at all K-points:
1295 : !! <k,n'|H(0)|k,n'> (hartree).
1296 : !! eigenq(nkpt_rbz*mband*nsppol) = 0-order eigenvalues at all shifted K-points:
1297 : !! <k+Q,n'|H(0)|k+Q,n'> (hartree).
1298 : !! eigen1(nkpt_rbz*2*nsppol*mband**2,3,mpert) = matrix of first-order:
1299 : !! <k+Q,n'|H(1)|k,n> (hartree) (calculated in dfpt_cgwf).
1300 : !! eig2nkq(2,mband*nsppol,nkpt,3,natom,3,natom*dim_eig2nkq) = second derivatives of
1301 : !! the electronic eigenvalues.
1302 : !! elph2_imagden = imaginary part of the denominator of the sum-over-state expression
1303 : !! for the electronic eigenenergy shift due to second-order electron-phonon
1304 : !! interation.
1305 : !! ieig2rf = integer for calculation type.
1306 : !! indsym(4,nsym,natom) = indirect indexing array for atom labels
1307 : !! (not used yet, but will be used with symmetries).
1308 : !! mband = maximum number of bands.
1309 : !! mpert = maximum number of perturbations.
1310 : !! natom = number of atoms in the unit cell.
1311 : !! npert = number of phonon perturbations, without taking into account directions:
1312 : !! natom.
1313 : !! nsym = number of symmetries (not used yet).
1314 : !! mpi_enreg = information about MPI parallelization.
1315 : !! nkpt_rbz = number of k-points for each perturbation.
1316 : !! nsppol = 1 for unpolarized, 2 for spin-polarized.
1317 : !! smdelta = integer controling the calculation of electron lifetimes.
1318 : !! symq(4,2,nsym) = 1 if symmetry preserves present qpoint. From littlegroup_q (not used yet).
1319 : !! symrec(3,3,nsym) = 3x3 matrices of the group symmetries (reciprocal space)
1320 : !! (not used yet).
1321 : !! symrel(3,3,nsym) = array containing the symmetries in real space (not used yet).
1322 : !! timrev = 1 if time-reversal preserves the q wavevector; 0 otherwise
1323 : !! (not in use yet).
1324 : !! dtset = OPTIONAL, dataset structure containing the input variable of the
1325 : !! calculation. This is required to use the k-interpolation routine.
1326 : !! eigenq_fine(mband_fine,mkpt_fine,nsppol_fine) = OPTIONAL, 0-order eigenvalues
1327 : !! at all shifted K-points: <k+Q,n'|H(0)|k+Q,n'> (hartree) of the
1328 : !! fine grid. This information is read from the WF dense k-grid file.
1329 : !! hdr_fine = OPTIONAL, header of the WF file of the fine k-point grid. This
1330 : !! variable is required for the k-interpolation routine.
1331 : !! hdr0 = header of the GS WF file of the corse k-point grid.
1332 : !!
1333 : !!
1334 : !! OUTPUT
1335 : !! eig2nkq(2,mband*nsppol,nkpt_rbz,3,npert,3,npert)= diagonal part of the
1336 : !! second-order eigenvalues: E^{(2),diag}_{k,q,j}.
1337 : !! eigbrd(2,mband*nsppol,nkpt_rbz,3,npert,3,npert)= OPTIONAL, array containing the
1338 : !! the contribution of each perturbations pair
1339 : !! to the eigenstate broadening (inverse lifetime)
1340 : !! computed statically (without phonon frequencies).
1341 : !!
1342 : !! SOURCE
1343 :
1344 108 : subroutine eig2tot(dtfil,xred,psps,pawtab,natom,bdeigrf,clflg,dim_eig2nkq,eigen0,eigenq,eigen1,eig2nkq,&
1345 54 : & elph2_imagden,esmear,ieig2rf,mband,mpert,npert,mpi_enreg,doccde,&
1346 54 : & nkpt_rbz,nsppol,smdelta,rprimd,dtset,occ_rbz,hdr0,eigbrd,eigenq_fine,hdr_fine)
1347 :
1348 : !Arguments ------------------------------------
1349 : !scalars
1350 : integer,intent(in) :: bdeigrf,dim_eig2nkq,ieig2rf,mband,mpert,natom,nkpt_rbz
1351 : integer,intent(in) :: npert,nsppol,smdelta
1352 : real(dp),intent(in) :: elph2_imagden,esmear
1353 : type(MPI_type),intent(inout) :: mpi_enreg
1354 : type(datafiles_type), intent(in) :: dtfil
1355 : type(pseudopotential_type), intent(inout) :: psps
1356 : !arrays
1357 : type(dataset_type), intent(in) :: dtset
1358 : integer,intent(in) :: clflg(3,mpert)
1359 : real(dp),intent(in) :: doccde(dtset%mband*dtset%nkpt*dtset%nsppol)
1360 : real(dp),intent(in) :: eigen0(nkpt_rbz*mband*nsppol)
1361 : real(dp),intent(in) :: eigen1(nkpt_rbz*2*nsppol*mband**2,3,mpert)
1362 : real(dp),intent(in) :: eigenq(nkpt_rbz*mband*nsppol)
1363 : real(dp),intent(in) :: occ_rbz(mband*nkpt_rbz*nsppol)
1364 : real(dp),intent(inout) :: eig2nkq(2,mband*nsppol,nkpt_rbz,3,npert,3,npert*dim_eig2nkq)
1365 : real(dp),intent(in) :: rprimd(3,3),xred(3,natom)
1366 : real(dp),intent(inout),optional :: eigbrd(2,mband*nsppol,nkpt_rbz,3,npert,3,npert)
1367 : real(dp),intent(in),pointer,optional :: eigenq_fine(:,:,:)
1368 : type(pawtab_type), intent(inout) :: pawtab(psps%ntypat*psps%usepaw)
1369 : type(hdr_type),intent(in) :: hdr0
1370 : type(hdr_type),intent(in),optional :: hdr_fine
1371 :
1372 : !Local variables-------------------------------
1373 : !tolerance for non degenerated levels
1374 : !scalars
1375 : integer :: band2tot_index,band_index,bantot,bandtot_index,iband,idir1,idir2
1376 : integer :: ikpt,ipert1,ipert2,isppol,jband,nkpt_sub,ikpt2,ncid
1377 : !integer :: ipw
1378 : character(len=fnlen) :: dscrpt,fname
1379 : integer :: master,me,spaceworld,ierr
1380 : integer :: mband_mem, mpert_
1381 : ! real(dp),parameter :: etol=1.0d-6
1382 : real(dp),parameter :: etol=1.0d-7
1383 : !real(dp),parameter :: etol=zero
1384 : real(dp) :: ar,ai,deltae,den,eig1_i1,eig1_i2,eigen_corr
1385 : real(dp) :: eig1_r1,eig1_r2,eig2_diai,den_av
1386 : real(dp) :: eig2_diar,eigbrd_i,eigbrd_r,wgt_int
1387 : !character(len=500) :: message
1388 : logical :: remove_inv,test_do_band
1389 2808 : type(crystal_t) :: Crystal
1390 54 : type(ebands_t) :: Bands
1391 : !type(eigr2d_t) :: eigr2d,eigi2d
1392 54 : type(fan_t) :: fan2d
1393 54 : type(gkk_t) :: gkk2d
1394 2808 : type(ddb_hdr_type) :: ddb_hdr
1395 54 : type(ddb_type) :: ddb
1396 : !arrays
1397 54 : integer,allocatable :: flg(:,:,:,:)
1398 54 : integer,allocatable :: nband_rbz(:)
1399 54 : integer,pointer :: kpt_fine_sub(:)
1400 : real(dp) :: tsec(2)
1401 54 : real(dp),allocatable :: center(:)
1402 108 : real(dp) :: eigen(mband*nsppol),eigen_prime(mband*nsppol)
1403 54 : real(dp),allocatable :: fan(:,:,:,:,:,:,:)
1404 54 : real(dp),allocatable :: gkk(:,:,:,:,:)
1405 54 : real(dp),allocatable :: smdfun(:,:)
1406 54 : real(dp),pointer :: wgt_sub(:)
1407 :
1408 : ! *********************************************************************
1409 :
1410 : !Init parallelism
1411 54 : master =0
1412 54 : spaceworld=mpi_enreg%comm_cell
1413 54 : me=mpi_enreg%me_kpt
1414 :
1415 : !Init interpolation method
1416 54 : if(present(eigenq_fine))then
1417 0 : ABI_MALLOC(center,(3))
1418 : end if
1419 :
1420 54 : call timab(148,1,tsec)
1421 :
1422 54 : if(nsppol==2)then
1423 1 : ABI_COMMENT('nsppol=2 is still under development. Be careful when using it ...')
1424 : end if
1425 :
1426 54 : band2tot_index =0
1427 54 : bandtot_index=0
1428 54 : band_index=0
1429 :
1430 : if(xmpi_paral==1) then
1431 270 : ABI_MALLOC(mpi_enreg%proc_distrb,(nkpt_rbz,mband,nsppol))
1432 162 : ABI_MALLOC(nband_rbz,(nkpt_rbz*nsppol))
1433 54 : ABI_SFREE(mpi_enreg%my_kpttab)
1434 162 : ABI_MALLOC(mpi_enreg%my_kpttab,(nkpt_rbz))
1435 : ! Assume the number of bands is the same for all k points.
1436 1096 : nband_rbz(:)=mband
1437 54 : call distrb2(mband,mband_mem,nband_rbz,nkpt_rbz,mpi_enreg%nproc_cell,nsppol,mpi_enreg)
1438 : end if
1439 :
1440 54 : if(ieig2rf == 4 ) then
1441 250 : ABI_MALLOC_OR_DIE(fan,(2*mband*nsppol,dtset%nkpt,3,natom,3,natom*dim_eig2nkq,mband), ierr)
1442 1604377 : fan(:,:,:,:,:,:,:) = zero
1443 25 : if(present(eigbrd))then
1444 229809 : eigbrd(:,:,:,:,:,:,:)=zero
1445 : end if
1446 : eigen_corr = 0
1447 : end if
1448 :
1449 54 : if(ieig2rf == 5 ) then
1450 232 : ABI_MALLOC_OR_DIE(gkk,(2*mband*nsppol,dtset%nkpt,3,natom,mband), ierr)
1451 1504589 : gkk(:,:,:,:,:) = zero
1452 29 : if(present(eigbrd))then
1453 1117545 : eigbrd(:,:,:,:,:,:,:)=zero
1454 : end if
1455 : eigen_corr = 0
1456 : end if
1457 :
1458 109 : do isppol=1,nsppol
1459 1097 : do ikpt =1,nkpt_rbz
1460 :
1461 1042 : if(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,mband,isppol,me)) then
1462 96 : band2tot_index = band2tot_index + 2*mband**2
1463 96 : bandtot_index = bandtot_index + mband
1464 96 : cycle
1465 : end if
1466 :
1467 946 : if(present(eigenq_fine))then
1468 0 : write(std_out,*) 'Start of the energy denominator interpolation method.'
1469 0 : nkpt_sub = 0
1470 : ! center is the k+q point around which we will average the kpt_fine
1471 0 : center = hdr0%kptns(:,ikpt)+ dtset%qptn(:)
1472 :
1473 0 : call kptfine_av(center,dtset%qptrlatt,hdr_fine%kptns,hdr_fine%nkpt,kpt_fine_sub,nkpt_sub,wgt_sub)
1474 : write(std_out,'(a,3f8.4,a,i3)') 'Number of k-points of the fine grid &
1475 0 : & around the k+Q point ',center,' is:',nkpt_sub
1476 0 : write(std_out,'(a,f10.5)') 'The sum of the weights of the k-points is: ',SUM(wgt_sub)
1477 : end if
1478 :
1479 946 : if(smdelta >0) then !broadening
1480 946 : if(.not.allocated(smdfun)) then
1481 1108 : ABI_MALLOC(smdfun,(mband,mband))
1482 : end if
1483 143646 : smdfun(:,:) = zero
1484 12086 : do iband=1,mband
1485 11140 : eigen(iband) = eigen0(iband+bandtot_index)
1486 12086 : eigen_prime(iband) =eigenq(iband+bandtot_index)
1487 : end do
1488 946 : if(esmear>tol6) then
1489 946 : call smeared_delta(eigen,eigen_prime,esmear,mband,smdelta,smdfun)
1490 : end if
1491 : end if
1492 :
1493 946 : ipert1=1 ! Suppose all perturbations lead to the same number of planewaves
1494 :
1495 12086 : do iband=1,bdeigrf
1496 :
1497 : ! If the k point and band belong to me, compute the contribution
1498 11140 : test_do_band=.true.
1499 11140 : if(mpi_enreg%proc_distrb(ikpt,iband,isppol)/=me)test_do_band=.false.
1500 :
1501 946 : if(test_do_band)then
1502 : ! ------------------------------------------------------------------------------------------------------!
1503 : ! ------- ieig2rf ==3 : Non dynamic traditional AHC theory with Sternheimer (computed in eig2stern.F90)-!
1504 : ! ------------------------------------------------------------------------------------------------------!
1505 : ! Note that ieig2rf==4 and ieig2rf==5 also goes into that part only for later printing of the ZPR in the ouput of abinit
1506 : ! later in the code
1507 11140 : if(ieig2rf==3 .or. ieig2rf==4 .or. ieig2rf==5) then
1508 33420 : do ipert1=1,npert
1509 100260 : do idir1=1,3
1510 66840 : if(clflg(idir1,ipert1)==0) cycle
1511 222800 : do ipert2=1,npert
1512 601560 : do idir2=1,3
1513 401040 : if(clflg(idir2,ipert2)==0)cycle
1514 : eig2_diar = zero ; eig2_diai = zero ; eigbrd_r = zero ; eigbrd_i = zero
1515 5137200 : do jband=1,mband
1516 4736160 : eig1_r1 = eigen1(2*jband-1+(iband-1)*2*mband+band2tot_index,idir1,ipert1)
1517 4736160 : eig1_r2 = eigen1(2*jband-1+(iband-1)*2*mband+band2tot_index,idir2,ipert2)
1518 4736160 : eig1_i1 = eigen1(2*jband+(iband-1)*2*mband+band2tot_index,idir1,ipert1)
1519 4736160 : eig1_i2 = - eigen1(2*jband+(iband-1)*2*mband+band2tot_index,idir2,ipert2) !the negative sign is from the CC
1520 : ! If no interpolation, fallback on to the previous
1521 : ! implementation
1522 4736160 : if(.not. present(eigenq_fine))then
1523 4736160 : deltae=eigenq(jband+bandtot_index)-eigen0(iband+bandtot_index)
1524 : end if
1525 4736160 : ar=eig1_r1*eig1_r2-eig1_i1*eig1_i2
1526 4736160 : ai=eig1_r1*eig1_i2+eig1_i1*eig1_r2
1527 :
1528 : ! Sum over all active space to retrieve the diagonal gauge
1529 : ! if(abs(deltae)>etol) then ! This is commented because
1530 : ! there is no problem with divergencies with elph2_imag != 0
1531 4736160 : if( present(eigenq_fine))then
1532 0 : den_av = zero
1533 0 : wgt_int = zero
1534 0 : do ikpt2=1,nkpt_sub
1535 : deltae=eigenq_fine(jband,kpt_fine_sub(ikpt2),1)&
1536 0 : & -eigen0(iband+bandtot_index)
1537 0 : den_av = den_av-(wgt_sub(ikpt2)*deltae)/(deltae**2+elph2_imagden**2)
1538 0 : wgt_int = wgt_int+wgt_sub(ikpt2)
1539 : end do
1540 0 : den = den_av/wgt_int
1541 : else
1542 4736160 : if(abs(elph2_imagden) < etol) then
1543 0 : if(abs(deltae)>etol) then
1544 0 : den=-one/(deltae**2+elph2_imagden**2)
1545 : else
1546 : den= zero
1547 : end if
1548 : else
1549 4736160 : den=-one/(deltae**2+elph2_imagden**2)
1550 : end if
1551 : end if
1552 :
1553 4736160 : if( present(eigenq_fine))then
1554 0 : eig2_diar=eig2_diar+ar*den
1555 0 : eig2_diai=eig2_diai+ai*den
1556 : else
1557 4736160 : eig2_diar=eig2_diar+ar*deltae*den
1558 4736160 : eig2_diai=eig2_diai+ai*deltae*den
1559 : end if
1560 :
1561 5137200 : if(present(eigbrd))then
1562 4736160 : if(smdelta >0) then !broadening
1563 4736160 : eigbrd_r = eigbrd_r + ar*smdfun(iband,jband)
1564 4736160 : eigbrd_i = eigbrd_i + ai*smdfun(iband,jband)
1565 : end if
1566 : end if
1567 : end do !jband
1568 :
1569 : ! Store the contribution
1570 : eig2nkq(1,iband+band_index,ikpt,idir1,ipert1,idir2,ipert2) = &
1571 401040 : & eig2nkq(1,iband+band_index,ikpt,idir1,ipert1,idir2,ipert2) + eig2_diar
1572 : eig2nkq(2,iband+band_index,ikpt,idir1,ipert1,idir2,ipert2) = &
1573 401040 : & eig2nkq(2,iband+band_index,ikpt,idir1,ipert1,idir2,ipert2) + eig2_diai
1574 :
1575 534720 : if(present(eigbrd))then
1576 401040 : if(smdelta >0) then !broadening
1577 401040 : eigbrd(1,iband+band_index,ikpt,idir1,ipert1,idir2,ipert2) = eigbrd_r
1578 401040 : eigbrd(2,iband+band_index,ikpt,idir1,ipert1,idir2,ipert2) = eigbrd_i
1579 : end if
1580 : end if
1581 :
1582 : end do !idir2
1583 : end do !ipert2
1584 : end do !idir1
1585 : end do !ipert1
1586 : end if !ieig2rf 3
1587 :
1588 : ! -------------------------------------------------------------------------------------------!
1589 : ! ------- ieig2rf ==4 Dynamic AHC using second quantization and Sternheimer from eig2stern -!
1590 : ! -------------------------------------------------------------------------------------------!
1591 11140 : if(ieig2rf ==4 ) then
1592 3264 : do ipert1=1,npert
1593 9792 : do idir1=1,3
1594 6528 : if(clflg(idir1,ipert1)==0) cycle
1595 21760 : do ipert2=1,npert
1596 58752 : do idir2=1,3
1597 39168 : if(clflg(idir2,ipert2)==0)cycle
1598 464640 : do jband=1,mband
1599 412416 : eig1_r1 = eigen1(2*jband-1+(iband-1)*2*mband+band2tot_index,idir1,ipert1)
1600 412416 : eig1_r2 = eigen1(2*jband-1+(iband-1)*2*mband+band2tot_index,idir2,ipert2)
1601 412416 : eig1_i1 = eigen1(2*jband+(iband-1)*2*mband+band2tot_index,idir1,ipert1)
1602 412416 : eig1_i2 = - eigen1(2*jband+(iband-1)*2*mband+band2tot_index,idir2,ipert2) !the negative sign is from the CC
1603 412416 : ar=eig1_r1*eig1_r2-eig1_i1*eig1_i2
1604 412416 : ai=eig1_r1*eig1_i2+eig1_i1*eig1_r2
1605 : ! Store the contribution
1606 : fan(2*iband-1+2*band_index,ikpt,idir1,ipert1,idir2,ipert2,jband) = &
1607 412416 : & fan(2*iband-1+2*band_index,ikpt,idir1,ipert1,idir2,ipert2,jband) + ar
1608 : fan(2*iband+2*band_index,ikpt,idir1,ipert1,idir2,ipert2,jband) = &
1609 451584 : & fan(2*iband+2*band_index,ikpt,idir1,ipert1,idir2,ipert2,jband) + ai
1610 : end do !jband
1611 : end do !idir2
1612 : end do !ipert2
1613 : end do !idir1
1614 : end do !ipert1
1615 : end if !ieig2rf 4
1616 : ! --------------------------------------------------------------------------------!
1617 : ! ------- ieig2rf ==5 Dynamic AHC with Sternheimer from eig2stern but print GKK -!
1618 : ! --------------------------------------------------------------------------------!
1619 11140 : if(ieig2rf ==5 ) then
1620 30156 : do ipert1=1,npert
1621 90468 : do idir1=1,3
1622 60312 : if(clflg(idir1,ipert1)==0) cycle
1623 801040 : do jband=1,mband
1624 720624 : eig1_r1 = eigen1(2*jband-1+(iband-1)*2*mband+band2tot_index,idir1,ipert1)
1625 720624 : eig1_i1 = eigen1(2*jband+(iband-1)*2*mband+band2tot_index,idir1,ipert1)
1626 : ! Store the contribution
1627 : gkk(2*iband-1+2*band_index,ikpt,idir1,ipert1,jband) = &
1628 720624 : & gkk(2*iband-1+2*band_index,ikpt,idir1,ipert1,jband) + eig1_r1
1629 : gkk(2*iband+2*band_index,ikpt,idir1,ipert1,jband) = &
1630 780936 : & gkk(2*iband+2*band_index,ikpt,idir1,ipert1,jband) + eig1_i1
1631 : end do !jband
1632 : end do !idir1
1633 : end do !ipert1
1634 : end if !ieig2rf 5
1635 : end if ! Selection of processor
1636 : end do !iband
1637 :
1638 946 : band2tot_index = band2tot_index + 2*mband**2
1639 946 : bandtot_index = bandtot_index + mband
1640 :
1641 1001 : if(present(eigenq_fine))then
1642 0 : ABI_FREE(kpt_fine_sub) ! Deallocate the variable
1643 0 : ABI_FREE(wgt_sub)
1644 : end if
1645 :
1646 : end do !ikpt
1647 109 : band_index = band_index + mband
1648 : end do !isppol
1649 :
1650 : !Accumulate eig2nkq and/or eigbrd
1651 : if(xmpi_paral==1) then
1652 54 : if(ieig2rf == 3) then
1653 0 : call xmpi_sum(eig2nkq,spaceworld,ierr)
1654 : end if
1655 54 : if(ieig2rf == 4) then
1656 25 : call xmpi_sum(eig2nkq,spaceworld,ierr)
1657 25 : call xmpi_sum(fan,spaceworld,ierr)
1658 : end if
1659 54 : if(ieig2rf == 5) then
1660 29 : call xmpi_sum(eig2nkq,spaceworld,ierr)
1661 29 : call xmpi_sum(gkk,spaceworld,ierr)
1662 : end if
1663 54 : if(present(eigbrd) .and. (ieig2rf == 3 .or. ieig2rf == 4 .or. ieig2rf == 5))then
1664 54 : if(smdelta >0) then
1665 54 : call xmpi_sum(eigbrd,spaceworld,ierr)
1666 : end if
1667 : end if
1668 54 : ABI_FREE(nband_rbz)
1669 54 : ABI_FREE(mpi_enreg%proc_distrb)
1670 54 : ABI_FREE(mpi_enreg%my_kpttab)
1671 : end if
1672 :
1673 54 : if (ieig2rf > 2) then
1674 : ! SP: Important: the second-order derivatives of the electronic energy printed here is the total one, including
1675 : ! both the Sternheimer and active-space part. Instead, the file _EIGR2D produced when ieig2rf > 2
1676 : ! contains ONLY the Sternheimer part and is therefore not the same values as here.
1677 : ! Note that in the case of ieig2rf == 1 or 2, the _EIGR2D file does contain the full one.
1678 54 : write(ab_out,'(a)')' Components of second-order derivatives of the electronic energy, EIGR2D.'
1679 54 : write(ab_out,'(a)')' For automatic tests, printing the matrix for the first k-point, first band, first atom.'
1680 54 : band_index = 0
1681 109 : do isppol=1,dtset%nsppol
1682 220 : do idir1=1,3
1683 715 : do idir2=1,3
1684 495 : ar=eig2nkq(1,1+band_index,1,idir1,1,idir2,1) ; if(abs(ar)<tol10)ar=zero
1685 495 : ai=eig2nkq(2,1+band_index,1,idir1,1,idir2,1) ; if(abs(ai)<tol10)ai=zero
1686 660 : write (ab_out,'(4i4,2es20.10)') idir1,1,idir2,1,ar,ai
1687 : end do ! idir2
1688 : end do ! idir1
1689 55 : band_index = band_index + mband
1690 109 : write(ab_out,'(a)')' '
1691 : end do
1692 : end if
1693 :
1694 54 : if(present(eigbrd))then
1695 54 : if(smdelta >0) then !broadening
1696 54 : write(ab_out,'(a)')' Components of second-order derivatives of the electronic energy, EIGI2D.'
1697 54 : write(ab_out,'(a)')' For automatic tests, printing the matrix for the first k-point, first band, first atom.'
1698 54 : band_index = 0
1699 109 : do isppol=1,dtset%nsppol
1700 220 : do idir1=1,3
1701 715 : do idir2=1,3
1702 495 : ar=eigbrd(1,1+band_index,1,idir1,1,idir2,1) ; if(abs(ar)<tol10)ar=zero
1703 495 : ai=eigbrd(2,1+band_index,1,idir1,1,idir2,1) ; if(abs(ai)<tol10)ai=zero
1704 660 : write (ab_out,'(4i4,2es20.10)') idir1,1,idir2,1,ar,ai
1705 : end do
1706 : end do
1707 55 : band_index = band_index + mband
1708 109 : write(ab_out,'(a)')' '
1709 : end do
1710 : end if
1711 : end if
1712 :
1713 54 : if(allocated(smdfun)) then
1714 54 : ABI_FREE(smdfun)
1715 : end if
1716 54 : if(present(eigenq_fine))then
1717 0 : ABI_FREE(center)
1718 : end if
1719 :
1720 : ! print _FAN file for this perturbation. Note that the Fan file will only be produced if
1721 : ! abinit is compiled with netcdf.
1722 :
1723 : ! Initialize crystal structure for FAN.nc and GKK.nc files
1724 54 : remove_inv=.false.
1725 54 : if(dtset%nspden==4 .and. dtset%usedmft==1) remove_inv=.true.
1726 : call crystal%init(dtset%amu_orig(:,1),dtset%spgroup,dtset%natom,dtset%npsp,psps%ntypat, &
1727 : & dtset%nsym,rprimd,dtset%typat,xred,dtset%ziontypat,dtset%znucl,1,&
1728 : & dtset%nspden==2.and.dtset%nsppol==1,remove_inv,hdr0%title,&
1729 108 : & dtset%symrel,dtset%tnons,dtset%symafm)
1730 : ! Electronic band energies.
1731 54 : bantot= dtset%mband*dtset%nkpt*dtset%nsppol
1732 : call bands%init(bantot, dtset%nelect,dtset%ne_qFD,dtset%nh_qFD,dtset%ivalence,&
1733 : & doccde,eigen0,hdr0%istwfk,hdr0%kptns,&
1734 : & hdr0%nband, hdr0%nkpt,hdr0%npwarr,hdr0%nsppol,hdr0%nspinor,&
1735 : & hdr0%tphysel,hdr0%tsmear,hdr0%occopt,hdr0%occ,hdr0%wtk,&
1736 : & hdr0%cellcharge, hdr0%kptopt, hdr0%kptrlatt_orig, hdr0%nshiftk_orig, hdr0%shiftk_orig, &
1737 54 : & hdr0%kptrlatt, hdr0%nshiftk, hdr0%shiftk)
1738 : !
1739 54 : if(ieig2rf == 4 ) then
1740 : ! Output of the Fan.nc file.
1741 25 : fname = strcat(dtfil%filnam_ds(4),"_FAN.nc")
1742 25 : call fan_init(fan,fan2d,dtset%mband,hdr0%nsppol,nkpt_rbz,dtset%natom)
1743 25 : NCF_CHECK_MSG(nctk_open_create(ncid, fname, xmpi_comm_self), "Creating FAN file")
1744 25 : NCF_CHECK(crystal%ncwrite(ncid))
1745 25 : NCF_CHECK(Bands%ncwrite(ncid))
1746 25 : call fan_ncwrite(fan2d,dtset%qptn(:),dtset%wtq, ncid)
1747 25 : NCF_CHECK(nf90_close(ncid))
1748 25 : ABI_FREE(fan)
1749 : end if
1750 : ! print _GKK.nc file for this perturbation. Note that the GKK file will only be produced if
1751 : ! abinit is compiled with netcdf.
1752 54 : if(ieig2rf == 5 ) then
1753 : ! Output of the GKK.nc file.
1754 29 : fname = strcat(dtfil%filnam_ds(4),"_GKK.nc")
1755 29 : call gkk_init(gkk,gkk2d,dtset%mband,hdr0%nsppol,nkpt_rbz,dtset%natom,3)
1756 29 : NCF_CHECK_MSG(nctk_open_create(ncid, fname, xmpi_comm_self), "Creating GKK file")
1757 29 : NCF_CHECK(crystal%ncwrite(ncid))
1758 29 : NCF_CHECK(bands%ncwrite(ncid))
1759 29 : call gkk_ncwrite(gkk2d,dtset%qptn(:),dtset%wtq, ncid)
1760 29 : NCF_CHECK(nf90_close(ncid))
1761 29 : ABI_FREE(gkk)
1762 : end if
1763 :
1764 : ! print _EIGI2D file for this perturbation
1765 54 : if (ieig2rf /= 5 ) then
1766 25 : if(smdelta>0) then
1767 :
1768 25 : mpert_ = dtset%natom
1769 125 : ABI_MALLOC(flg,(3,mpert_,3,mpert_))
1770 1425 : flg = one
1771 : call ddb%init(dtset, 1, mpert_, &
1772 : mband=bdeigrf,&
1773 : nkpt=nkpt_rbz,&
1774 : kpt=dtset%kptns(:,1:nkpt_rbz),&
1775 25 : with_d2eig=.true.)
1776 :
1777 : ! Create the ddb header
1778 25 : dscrpt=' Note : temporary (transfer) database '
1779 : call ddb_hdr%init(dtset,psps,pawtab,dscrpt,1,&
1780 : mpert=mpert_,&
1781 : xred=xred,occ=occ_rbz,&
1782 : mband=bdeigrf / dtset%nsppol,&
1783 : nkpt=nkpt_rbz,&
1784 25 : kpt=dtset%kptns(:,1:nkpt_rbz))
1785 :
1786 25 : call ddb%set_qpt(1, dtset%qptn)
1787 25 : call ddb%set_d2eig_reshape(1, eigbrd, flg, blktyp=BLKTYP_d2eig_im)
1788 :
1789 25 : call ddb_hdr%set_typ(ddb%nblok, ddb%typ)
1790 :
1791 25 : call ddb_hdr%open_write(dtfil%fnameabo_eigi2d, with_psps=1)
1792 25 : call ddb%write_d2eig(ddb_hdr, 1)
1793 :
1794 25 : call ddb_hdr%close()
1795 25 : call ddb_hdr%free()
1796 25 : call ddb%free()
1797 :
1798 25 : ABI_FREE(flg)
1799 :
1800 : end if !smdelta
1801 : end if
1802 : !end if ! master
1803 :
1804 54 : ABI_SFREE(fan)
1805 54 : ABI_SFREE(gkk)
1806 :
1807 54 : call crystal%free()
1808 54 : call Bands%free()
1809 54 : call fan_free(fan2d)
1810 54 : call gkk_free(gkk2d)
1811 :
1812 54 : call timab(148,2,tsec)
1813 :
1814 162 : end subroutine eig2tot
1815 : !!***
1816 :
1817 : !!****f* m_eig2d/outbsd
1818 : !! NAME
1819 : !! outbsd
1820 : !!
1821 : !! FUNCTION
1822 : !! output bsd file for one perturbation (used for elphon calculations in anaddb)
1823 : !!
1824 : !! INPUTS
1825 : !! bdeigrf=number of bands for which the derivatives of the eigenvalues have been computed
1826 : !! dtset = dataset variable for run flags
1827 : !! eig2nkq= second ordre eigenvalue (or electron lifetime) that must be printed out
1828 : !! mpert= maximum number of perturbations
1829 : !! nkpt_rbz= number of k-points for perturbation
1830 : !! unitout= writting unit of file
1831 : !!
1832 : !! OUTPUTS
1833 : !! to file
1834 : !!
1835 : !! NOTE
1836 : !! This function is deprecated. One should write through ddb object instead.
1837 : !!
1838 : !! SOURCE
1839 :
1840 0 : subroutine outbsd(bdeigrf,dtset,eig2nkq,mpert,nkpt_rbz,unitout)
1841 :
1842 : !Arguments ------------------------------------
1843 : !scalars
1844 : integer,intent(in) :: bdeigrf,mpert,nkpt_rbz,unitout
1845 : type(dataset_type),intent(in) :: dtset
1846 : !arrays
1847 : real(dp),intent(in) :: eig2nkq(2,dtset%mband*dtset%nsppol,nkpt_rbz,3,mpert,3,mpert)
1848 :
1849 : !Local variables -------------------------
1850 : !scalars
1851 : integer :: bandtot_index,iband,idir1,idir2,ikpt,ipert1,ipert2,isppol
1852 :
1853 : ! *********************************************************************
1854 :
1855 :
1856 : !output information in this file
1857 0 : write(unitout,*)
1858 0 : write(unitout,'(a,i8)') ' 2nd eigenvalue derivatives - # elements :', 9*dtset%natom**2
1859 0 : write(unitout,'(a,3es16.8,a)') ' qpt', dtset%qptn(:), ' 1.0'
1860 :
1861 : !output RF eigenvalues
1862 :
1863 0 : do ikpt=1,nkpt_rbz
1864 : ! bandtot_index differs from zero only in the spin-polarized case
1865 0 : bandtot_index=0
1866 0 : write (unitout,'(a,3es16.8)') ' K-point:', dtset%kptns(:,ikpt)
1867 0 : do isppol=1,dtset%nsppol
1868 0 : do iband=1,bdeigrf
1869 0 : write (unitout,'(a,i5)') ' Band:', iband+bandtot_index
1870 : ! write (unitout,*) 'ipert1 ','idir1 ','ipert2 ','idir2 ','Real ','Im '
1871 0 : do ipert2=1,mpert
1872 0 : do idir2=1,3
1873 0 : do ipert1=1,mpert
1874 0 : do idir1=1,3
1875 0 : write (unitout,'(4i4,2d22.14)') idir1,ipert1,idir2,ipert2,&
1876 0 : & eig2nkq(1,iband+bandtot_index,ikpt,idir1,ipert1,idir2,ipert2),&
1877 0 : & eig2nkq(2,iband+bandtot_index,ikpt,idir1,ipert1,idir2,ipert2)
1878 : end do !idir2
1879 : end do !ipert2
1880 : end do !idir1
1881 : end do !ipert1
1882 : end do !iband
1883 0 : bandtot_index = bandtot_index + dtset%mband
1884 : end do !isppol
1885 : end do !ikpt
1886 :
1887 : !close bsd file
1888 0 : close (unitout)
1889 :
1890 0 : end subroutine outbsd
1891 : !!***
1892 :
1893 : !!****f* m-eig2d/smeared_delta
1894 : !! NAME
1895 : !! smeared_delta
1896 : !!
1897 : !! FUNCTION
1898 : !! This subroutine calculates the smeared delta that weights matrix elements:
1899 : !! \delta (\epsilon_{kn}-\epsilon_{k+Q,n'})
1900 : !!
1901 : !! INPUTS
1902 : !! eigen0(mband*nsppol) : Eigenvalues at point K
1903 : !! eigenq(mband*nsppol) : Eigenvalues at point K+Q
1904 : !! mband : maximum number of bands
1905 : !! smdelta : Variable controlling the smearinf scheme
1906 : !!
1907 : !! OUTPUT
1908 : !! smdfunc(mband,mband) : Smeared delta function weight corresponding to \delta(\epsilon_{n,k} - \epsilon_{n',k+Q})
1909 : !!
1910 : !! SOURCE
1911 :
1912 2143 : subroutine smeared_delta(eigen0,eigenq,esmear,mband,smdelta,smdfunc)
1913 :
1914 : !Arguments ------------------------------------
1915 : !scalars
1916 : integer,intent(in) :: mband,smdelta
1917 : !arrays
1918 : real(dp),intent(in) :: eigen0(mband),eigenq(mband),esmear
1919 : real(dp),intent(out) :: smdfunc(mband,mband)
1920 :
1921 : !Local variables-------------------------------
1922 : !tolerance for non degenerated levels
1923 : !scalars
1924 : integer :: ii,jj
1925 : real(dp) :: aa,dsqrpi,gauss,xx
1926 : character(len=500) :: message
1927 :
1928 : ! *********************************************************************
1929 :
1930 :
1931 : !---------------------------------------------------------
1932 : !Ordinary (unique) smearing function
1933 : !---------------------------------------------------------
1934 :
1935 2143 : if(smdelta==1)then
1936 :
1937 : ! Fermi-Dirac
1938 26119 : do ii=1,mband
1939 302719 : do jj= 1,mband
1940 276600 : xx= ( eigen0(ii) - eigenq(jj) )/esmear
1941 300576 : smdfunc(ii,jj)=0.25_dp/esmear/(cosh(xx/2.0_dp))**2
1942 : end do
1943 : end do
1944 :
1945 0 : else if(smdelta==2 .or. smdelta==3)then
1946 :
1947 : ! Cold smearing of Marzari, two values of the "a" parameter being possible
1948 : ! first value gives minimization of the bump
1949 : if(smdelta==2)aa=-.5634
1950 : ! second value gives monotonic occupation function
1951 0 : if(smdelta==3)aa=-.8165
1952 :
1953 0 : dsqrpi=1.0_dp/sqrt(pi)
1954 0 : do ii=1,mband
1955 0 : do jj=1,mband
1956 0 : xx= ( eigen0(ii) - eigenq(jj) ) / esmear
1957 0 : gauss=dsqrpi*exp(-xx**2)/esmear
1958 0 : smdfunc(ii,jj)=gauss*(1.5_dp+xx*(-aa*1.5_dp+xx*(-1.0_dp+aa*xx)))
1959 : end do
1960 : end do
1961 :
1962 0 : else if(smdelta==4)then
1963 :
1964 : ! First order Hermite-Gaussian of Paxton and Methfessel
1965 0 : dsqrpi=1.0_dp/sqrt(pi)
1966 0 : do ii=1,mband
1967 0 : do jj=1,mband
1968 0 : xx= ( eigen0(ii) - eigenq (jj) ) / esmear
1969 0 : smdfunc(ii,jj)=dsqrpi*(1.5_dp-xx**2)*exp(-xx**2)/esmear
1970 : end do
1971 : end do
1972 :
1973 0 : else if(smdelta==5)then
1974 :
1975 : ! Gaussian smearing
1976 0 : dsqrpi=1.0_dp/sqrt(pi)
1977 0 : do ii=1,mband
1978 0 : do jj=1,mband
1979 0 : xx= ( eigen0(ii) - eigenq (jj) ) / esmear
1980 0 : smdfunc(ii,jj)=dsqrpi*exp(-xx**2)/esmear
1981 : end do
1982 : end do
1983 :
1984 : else
1985 0 : write(message, '(a,i0,a)' )' Smdelta= ',smdelta,' is not allowed in smdfunc'
1986 0 : ABI_BUG(message)
1987 : end if
1988 :
1989 2143 : end subroutine smeared_delta
1990 : !!***
1991 :
1992 : !!****f* m_eig2d/elph2_fanddw
1993 : !! NAME
1994 : !! elph2_fanddw
1995 : !!
1996 : !! FUNCTION
1997 : !! This routine calculates the zero-point motion corrections
1998 : !! due to the Fan term or to the DDW term..
1999 : !!
2000 : !! INPUTS
2001 : !! dim_eig2nkq=1 if eig2nkq is to be computed
2002 : !! displ(2*3*natom*3*natom)=the displacements of atoms in cartesian coordinates.
2003 : !! eig2nkq(2,mband*nsppol,nkpt,3,natom,3,natom*dim_eig2nkq)=one half second derivatives of the electronic eigenvalues
2004 : !! gprimd(3,3)=dimensional primitive translations for reciprocal space ($\textrm{bohr}^{-1}$)
2005 : !! mband= maximum number of bands
2006 : !! natom= number of atoms in the unit cell
2007 : !! nkpt= number of k-points
2008 : !! nsppol= 1 for unpolarized, 2 for spin-polarized
2009 : !! option 1 for Fan term, 2 for DDW term
2010 : !! phfrq(3*natom)=phonon frequencies
2011 : !! (prtvol > 4) if the mode decomposition is to be printed
2012 : !!
2013 : !! OUTPUT
2014 : !! eigen_corr(mband*nkpt*nsppol)= T=0 correction to the electronic eigenvalues, due to the Fan term.
2015 : !!
2016 : !! SOURCE
2017 :
2018 97 : subroutine elph2_fanddw(dim_eig2nkq,displ,eig2nkq,eigen_corr,gprimd,mband,natom,nkpt,nsppol,option,phfrq,prtvol)
2019 :
2020 : !Arguments ------------------------------------
2021 : !scalars
2022 : integer,intent(in) :: dim_eig2nkq,mband,natom,nkpt,nsppol,option,prtvol
2023 :
2024 : !arrays
2025 : real(dp) :: gprimd(3,3)
2026 : real(dp),intent(in) :: displ(2*3*natom*3*natom)
2027 : real(dp),intent(in) :: eig2nkq(2,mband*nsppol,nkpt,3,natom,3,natom*dim_eig2nkq)
2028 : real(dp),intent(in) :: phfrq(3*natom)
2029 : real(dp),intent(out) :: eigen_corr(mband*nkpt*nsppol)
2030 :
2031 : !Local variables-------------------------------
2032 : !scalars
2033 : integer,parameter :: neigs_per_line=6
2034 : integer :: iatom1,iatom2,idir1,idir2,iband,ikpt,imode,index,isppol, imin, ii
2035 : real(dp) :: d_at1_dir1_re,d_at1_dir1_im
2036 : real(dp) :: d_at1_dir2_re,d_at1_dir2_im
2037 : real(dp) :: d_at2_dir1_re,d_at2_dir1_im
2038 : real(dp) :: d_at2_dir2_re,d_at2_dir2_im
2039 : real(dp) :: e2_im,e2_re
2040 97 : real(dp), allocatable :: eigen_corr_mode(:)
2041 : character(len=500) :: message
2042 : character(len=20) :: eig_format, line_format
2043 : !arrays
2044 : real(dp) :: displ2cart(2,3,3),displ2red(2,3,3),tmp_displ2(2,3,3)
2045 :
2046 : ! *********************************************************************
2047 :
2048 :
2049 97 : if(option/=1 .and. option/=2)then
2050 0 : write(message,'(a,i0)')' The argument option should be 1 or 2, while it is found that option=',option
2051 0 : ABI_BUG(message)
2052 : end if
2053 :
2054 : !printing options
2055 97 : eig_format='f16.8'
2056 97 : write(line_format,'(a,i1,a6,a)') '(',neigs_per_line,eig_format,')'
2057 :
2058 97 : if (prtvol > 4) then
2059 0 : write(message,'(a,a)')ch10,' ================================================================================'
2060 0 : call wrtout(ab_out_default,message,'COLL')
2061 0 : if (option==1) then
2062 0 : write(message,'(a)') ' ---- Begin Fan contributions to eigenvalues renormalization by mode ----'
2063 0 : call wrtout(ab_out_default,message,'COLL')
2064 0 : else if (option==2) then
2065 0 : write(message,'(a)') ' ---- Begin DDW contributions to eigenvalues renormalization by mode ----'
2066 0 : call wrtout(ab_out_default,message,'COLL')
2067 : end if
2068 : end if
2069 :
2070 291 : ABI_MALLOC(eigen_corr_mode,(mband*nkpt*nsppol))
2071 :
2072 16533 : eigen_corr(:)=zero
2073 667 : do imode=1,3*natom
2074 98802 : eigen_corr_mode(:)=zero
2075 :
2076 570 : if (phfrq(imode)>tol6) then
2077 1452 : do iatom1=1,natom
2078 5502 : do iatom2=1,natom
2079 :
2080 16200 : do idir1=1,3
2081 52650 : do idir2=1,3
2082 : ! Compute the mean cartesian displacements
2083 36450 : d_at1_dir1_re=displ(1 + 2*(idir1-1 +3*(iatom1-1 +natom*(imode-1))))
2084 36450 : d_at1_dir1_im=displ(2 + 2*(idir1-1 +3*(iatom1-1 +natom*(imode-1))))
2085 36450 : d_at2_dir2_re=displ(1 + 2*(idir2-1 +3*(iatom2-1 +natom*(imode-1))))
2086 36450 : d_at2_dir2_im=displ(2 + 2*(idir2-1 +3*(iatom2-1 +natom*(imode-1))))
2087 :
2088 48600 : if(option==1)then
2089 : ! Compute the mean displacement correlation at T=0.
2090 : ! Consistent with Eqs.(7) and (8) of PRB51, 8610 (1995) [[cite:Lee1995]], specialized for the contribution of one q point.
2091 : ! but generalized to two different atoms. Note that the complex conjugate is taken on the second direction.
2092 : displ2cart(1,idir1,idir2)=(d_at1_dir1_re*d_at2_dir2_re+ &
2093 22005 : & d_at1_dir1_im*d_at2_dir2_im )/(two*phfrq(imode))
2094 : displ2cart(2,idir1,idir2)=(d_at1_dir1_im*d_at2_dir2_re- &
2095 22005 : & d_at1_dir1_re*d_at2_dir2_im )/(two*phfrq(imode))
2096 14445 : else if(option==2)then
2097 : ! Compute the mean square displacement correlation of each atom at T=0, and take mean over iatom1 and iatom2.
2098 : ! See Eqs.(7) and (8) of PRB51, 8610 (1995) [[cite:Lee1995]], specialized for the contribution of one q point.
2099 : ! Note that the complex conjugate is taken on the second direction.
2100 : ! Also, note the overall negative sign, to make it opposite to the Fan term.
2101 14445 : d_at1_dir2_re=displ(1 + 2*(idir2-1 +3*(iatom1-1 +natom*(imode-1))))
2102 14445 : d_at1_dir2_im=displ(2 + 2*(idir2-1 +3*(iatom1-1 +natom*(imode-1))))
2103 14445 : d_at2_dir1_re=displ(1 + 2*(idir1-1 +3*(iatom2-1 +natom*(imode-1))))
2104 14445 : d_at2_dir1_im=displ(2 + 2*(idir1-1 +3*(iatom2-1 +natom*(imode-1))))
2105 : displ2cart(1,idir1,idir2)=-(d_at1_dir1_re*d_at1_dir2_re+ &
2106 : & d_at1_dir1_im*d_at1_dir2_im+ &
2107 : & d_at2_dir1_re*d_at2_dir2_re+ &
2108 14445 : & d_at2_dir1_im*d_at2_dir2_im )/(four*phfrq(imode))
2109 : displ2cart(2,idir1,idir2)=-(d_at1_dir1_im*d_at1_dir2_re- &
2110 : & d_at1_dir1_re*d_at1_dir2_im+ &
2111 : & d_at2_dir1_im*d_at2_dir2_re- &
2112 14445 : & d_at2_dir1_re*d_at2_dir2_im )/(four*phfrq(imode))
2113 : end if
2114 : end do
2115 : end do
2116 : ! Switch to reduced coordinates in two steps
2117 4050 : tmp_displ2(:,:,:)=zero
2118 16200 : do idir1=1,3
2119 52650 : do idir2=1,3
2120 376650 : tmp_displ2(:,:,idir1)=tmp_displ2(:,:,idir1)+displ2cart(:,:,idir2)*gprimd(idir2,idir1)
2121 : end do
2122 : end do
2123 4050 : displ2red(:,:,:)=zero
2124 16200 : do idir1=1,3
2125 52650 : do idir2=1,3
2126 376650 : displ2red(:,idir1,:)=displ2red(:,idir1,:)+tmp_displ2(:,idir2,:)*gprimd(idir2,idir1)
2127 : end do
2128 : end do
2129 : ! Compute the T=0 shift due to this q point
2130 17238 : do idir1=1,3
2131 52650 : do idir2=1,3
2132 296424 : do ikpt=1,nkpt
2133 532170 : do isppol=1,nsppol
2134 3559464 : do iband=1,mband
2135 3063744 : index=iband+mband*(isppol-1 + nsppol*(ikpt-1))
2136 3063744 : e2_re=eig2nkq(1,iband+mband*(isppol-1),ikpt,idir1,iatom1,idir2,iatom2)
2137 3063744 : e2_im=eig2nkq(2,iband+mband*(isppol-1),ikpt,idir1,iatom1,idir2,iatom2)
2138 : eigen_corr(index)=eigen_corr(index)+&
2139 3063744 : & e2_re*displ2red(1,idir1,idir2)-e2_im*displ2red(2,idir1,idir2)
2140 : eigen_corr_mode(index)=eigen_corr_mode(index)+&
2141 3311640 : & e2_re*displ2red(1,idir1,idir2)-e2_im*displ2red(2,idir1,idir2)
2142 : end do ! band
2143 : end do ! spin
2144 : end do ! kpt
2145 : end do ! dir2
2146 : end do ! dir1
2147 : end do ! atom2
2148 : end do ! atom1
2149 : end if
2150 :
2151 667 : if (prtvol > 4) then
2152 : ! Print the corrections by mode
2153 0 : write(message,'(a,i1)') ' imode= ',imode
2154 0 : call wrtout(ab_out_default,message,'COLL')
2155 :
2156 0 : do ikpt=1,nkpt
2157 0 : do isppol=1,nsppol
2158 0 : write(message,'(a,i4,a,i1)')' ikpt= ',ikpt,' ispin= ',isppol
2159 0 : call wrtout(ab_out_default,message,'COLL')
2160 :
2161 0 : imin = mband * (isppol-1 + nsppol*(ikpt-1))
2162 0 : do ii=0, (mband-1)/neigs_per_line
2163 0 : write(message, line_format) (eigen_corr_mode(iband+imin), &
2164 0 : & iband = 1 + ii * neigs_per_line, min(mband, (ii+1)*neigs_per_line))
2165 0 : call wrtout(ab_out_default,message,'COLL')
2166 : end do
2167 : end do
2168 : end do
2169 : end if
2170 :
2171 : end do ! mode
2172 :
2173 97 : if (prtvol > 4) then
2174 0 : if (option==1) then
2175 0 : write(message,'(a)') ' ---- End Fan contribution to eigenvalues renormalization by mode ----'
2176 0 : call wrtout(ab_out_default,message,'COLL')
2177 0 : else if (option==2) then
2178 0 : write(message,'(a)') ' ---- End DDW contribution to eigenvalues renormalization by mode ----'
2179 0 : call wrtout(ab_out_default,message,'COLL')
2180 : end if
2181 0 : write(message,'(a,a)')' ================================================================================', ch10
2182 0 : call wrtout(ab_out_default,message,'COLL')
2183 : end if
2184 :
2185 97 : ABI_FREE(eigen_corr_mode)
2186 :
2187 97 : end subroutine elph2_fanddw
2188 : !!***
2189 :
2190 0 : END MODULE m_eig2d
2191 : !!***
|
