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