Line data Source code
1 : !!****m* ABINIT/m_ephtk
2 : !! NAME
3 : !! m_ephtk
4 : !!
5 : !! FUNCTION
6 : !! Helper functions common to e-ph calculations.
7 : !!
8 : !! COPYRIGHT
9 : !! Copyright (C) 2008-2026 ABINIT group (MG)
10 : !! This file is distributed under the terms of the
11 : !! GNU General Public License, see ~abinit/COPYING
12 : !! or http://www.gnu.org/copyleft/gpl.txt .
13 : !!
14 : !! SOURCE
15 :
16 : #if defined HAVE_CONFIG_H
17 : #include "config.h"
18 : #endif
19 :
20 : #include "abi_common.h"
21 :
22 : module m_ephtk
23 :
24 : use defs_basis
25 : use m_abicore
26 : use m_errors
27 : use m_xmpi
28 : use m_dtset
29 : use m_ebands
30 : use m_crystal
31 : use m_krank
32 : use m_xmpi
33 :
34 : use m_fstrings, only : itoa, sjoin, ltoa, ftoa, ktoa
35 : use m_bz_mesh, only : isamek
36 : use m_fftcore, only : get_kg
37 :
38 : implicit none
39 :
40 : private
41 :
42 : public :: ephtk_set_phmodes_skip ! Setup a mask to skip accumulating the contribution of certain phonon modes.
43 : public :: ephtk_set_pertables ! Set tables for parallelism over perturbations from my_npert and comm
44 : public :: ephtk_mkqtabs ! Build tables with correspondence between q-points as needed by complete_gamma.
45 : public :: ephtk_gam_atm2qnu ! Compute phonon linewidths from gamma matrix in reduced coordinates.
46 : public :: ephtk_gkknu_from_atm ! Transform the gkk matrix elements from (atom, red_direction) basis to phonon-mode basis.
47 : public :: ephtk_update_ebands ! Update ebands according to dtset%occopt, tsmear, mbpt_sciss, eph_fermie, eph_extrael
48 : public :: ephtk_get_mpw_gmax ! Compute maximum number of plane-waves over k and k+q where k and k+q are in the BZ.
49 : public :: ephtk_v1atm_to_vqnu ! Receive potentials in atomic representation and return potential in phonon representation
50 : public :: ephtk_skip_phmode ! Ignore contribution of phonon mode depending on phonon frequency value or mode index.
51 : !!***
52 :
53 : real(dp),public,parameter :: EPHTK_WTOL = tol6
54 : ! Tolerance for phonon frequencies to be ignored.
55 : ! Lambda coefficients are set to zero when abs(w) < EPHTK_WTOL
56 : ! This tolerance is also used in the integrals of a2F(w).
57 :
58 : contains !=====================================================
59 : !!***
60 :
61 : !!****f* m_ephtk/ephtk_set_phmodes_skip
62 : !! NAME
63 : !! ephtk_set_phmodes_skip
64 : !!
65 : !! FUNCTION
66 : !! Setup a mask to skip accumulating the contribution of certain phonon modes.
67 : !!
68 : !! INPUT
69 : !! eph_phrange=Abinit input variable.
70 : !!
71 : !! OUTPUT
72 : !! phmodes_skip(natom3) For each mode: 1 to skip the contribution given by this phonon branch else 0
73 : !!
74 : !! SOURCE
75 :
76 65 : subroutine ephtk_set_phmodes_skip(natom, eph_phrange, phmodes_skip)
77 :
78 : !Arguments ------------------------------------
79 : integer,intent(in) :: natom
80 : !arrays
81 : integer,intent(in) :: eph_phrange(2)
82 : integer,allocatable,intent(out) :: phmodes_skip(:)
83 :
84 : !Local variables ------------------------------
85 : integer :: natom3
86 : ! *************************************************************************
87 :
88 : ! Setup a mask to skip accumulating the contribution of certain phonon modes.
89 : ! By default do not skip, if set skip all but specified
90 65 : natom3 = natom * 3
91 195 : ABI_MALLOC(phmodes_skip, (natom3))
92 452 : phmodes_skip = 0
93 :
94 67 : if (all(eph_phrange /= 0)) then
95 : if (minval(abs(eph_phrange)) < 1 .or. &
96 5 : maxval(abs(eph_phrange)) > natom3 .or. &
97 : abs(eph_phrange(2)) < abs(eph_phrange(1))) then
98 0 : ABI_ERROR('Invalid range for eph_phrange. Should be between [1, 3*natom] and eph_modes(2) > eph_modes(1)')
99 : end if
100 3 : if (all(eph_phrange > 0)) then
101 1 : call wrtout(std_out, sjoin(" Including phonon modes between [", itoa(eph_phrange(1)), ',', itoa(eph_phrange(2)), "]"))
102 7 : phmodes_skip = 1
103 7 : phmodes_skip(eph_phrange(1):eph_phrange(2)) = 0
104 0 : else if (all(eph_phrange < 0)) then
105 : call wrtout(std_out, sjoin(" Excluding phonon modes between [", &
106 0 : itoa(abs(eph_phrange(1))), ',', itoa(abs(eph_phrange(2))), "]"))
107 0 : phmodes_skip = 0
108 0 : phmodes_skip(abs(eph_phrange(1)):abs(eph_phrange(2))) = 1
109 : else
110 0 : ABI_ERROR(sjoin("Invalid eph_phrange: ", itoa(eph_phrange(1)), ',', itoa(eph_phrange(2))))
111 : end if
112 : end if
113 :
114 65 : end subroutine ephtk_set_phmodes_skip
115 : !!***
116 :
117 : !!****f* m_ephtk/ephtk_set_pertables
118 : !! NAME
119 : !! ephtk_set_pertables
120 : !!
121 : !! FUNCTION
122 : !! Build tables for parallelism over perturbations from my_npert and comm
123 : !!
124 : !! INPUT
125 : !! natom: Number of atoms
126 : !! my_npert: Number of atomic perturbations or phonon modes treated by this MPI rank.
127 : !! comm: MPI communicator for parallelism over atomic perturbations.
128 : !!
129 : !! OUTPUT
130 : !! integer,allocatable :: my_pinfo(:,:)
131 : !! my_pinfo(3, my_npert)
132 : !! my_pinfo(1, ip) gives the `idir` index of the ip-th perturbation.
133 : !! my_pinfo(2, ip) gives the `ipert` index of the ip-th perturbation.
134 : !! my_pinfo(3, ip) gives `pertcase`=idir + (ipert-1)*3
135 : !! integer,allocatable :: pert_table(:,:)
136 : !! pert_table(2, natom3)
137 : !! pert_table(1, npert): rank of the processor treating this atomic perturbation.
138 : !! pert_table(2, npert): imyp index in my_pinfo table, -1 if this rank is not treating ipert.
139 : !!
140 : !! SOURCE
141 :
142 58 : subroutine ephtk_set_pertables(natom, my_npert, pert_table, my_pinfo, comm)
143 :
144 : !Arguments ------------------------------------
145 : integer,intent(in) :: natom, my_npert, comm
146 : !arrays
147 : integer,allocatable,intent(out) :: pert_table(:,:), my_pinfo(:,:)
148 :
149 : !Local variables ------------------------------
150 : !scalars
151 : integer :: iatom, idir, pertcase, bstart, bstop, ii, ip, natom3, my_rank, nproc
152 : !arrays
153 116 : integer :: all_pinfo(3, natom*3)
154 : ! *************************************************************************
155 :
156 58 : my_rank = xmpi_comm_rank(comm); nproc = xmpi_comm_size(comm)
157 :
158 : ! Build table with list of perturbations treated by this CPU.
159 58 : natom3 = natom * 3
160 174 : ABI_MALLOC(my_pinfo, (3, my_npert))
161 174 : ABI_MALLOC(pert_table, (2, natom3))
162 :
163 173 : do iatom=1,natom
164 518 : do idir=1,3
165 345 : pertcase = idir + (iatom-1) * 3
166 1380 : all_pinfo(:, pertcase) = [idir, iatom, pertcase]
167 460 : pert_table(1, pertcase) = (pertcase - 1) / (natom3 / nproc)
168 : end do
169 : end do
170 58 : bstart = (natom3 / nproc) * my_rank + 1
171 58 : bstop = bstart + my_npert - 1
172 1496 : my_pinfo = all_pinfo(:, bstart:bstop)
173 :
174 403 : pert_table(2, :) = -1
175 403 : do ii=1,my_npert
176 345 : ip = my_pinfo(3, ii)
177 403 : pert_table(2, ip) = ii
178 : end do
179 : !write(std_out,*)"my_npert", my_npert, "nproc", nproc; write(std_out,*)"my_pinfo", my_pinfo
180 :
181 58 : end subroutine ephtk_set_pertables
182 : !!***
183 :
184 : !!****f* m_ephtk/ephtk_mkqtabs
185 : !! NAME
186 : !! ephtk_mkqtabs
187 : !!
188 : !! FUNCTION
189 : !! Build tables with correspondence between q-points in the IBZ/BZ as needed by complete_gamma.
190 : !!
191 : !! INPUT
192 : !! cryst<crystal_t>=Crystal structure.
193 : !! nqibz, qibz = Points in the IBZ
194 : !! nqbz, qbz = Points in the BZ
195 : !!
196 : !! OUTPUT
197 : !! qirredtofull(nqibz) = mapping irred to full qpoints
198 : !! qpttoqpt(2, cryst%nsym, nqbz)) = qpoint index mapping under symops.
199 : !!
200 : !! SOURCE
201 :
202 10 : subroutine ephtk_mkqtabs(cryst, nqibz, qibz, nqbz, qbz, qirredtofull, qpttoqpt)
203 :
204 : !Arguments ------------------------------------
205 : type(crystal_t),intent(in) :: cryst
206 : integer,intent(in) :: nqibz, nqbz
207 : !arrays
208 : real(dp),intent(in) :: qibz(3, nqibz), qbz(3, nqbz)
209 : integer,allocatable :: qirredtofull(:),qpttoqpt(:,:,:)
210 :
211 : !Local variables ------------------------------
212 : !scalars
213 : integer :: iq_bz, iq_ibz, isq_bz, isym
214 10 : type(krank_t) :: qrank
215 : !arrays
216 : integer :: g0(3)
217 : real(dp) :: qirr(3), tmp_qpt(3)
218 : ! *************************************************************************
219 :
220 10 : call qrank%init(nqbz, qbz)
221 :
222 : ! Compute index of IBZ q-point in the BZ array
223 152 : ABI_CALLOC(qirredtofull, (nqibz))
224 :
225 132 : do iq_ibz=1,nqibz
226 488 : qirr = qibz(:,iq_ibz)
227 122 : iq_bz = qrank%get_index(qirr)
228 132 : if (iq_bz /= -1) then
229 122 : ABI_CHECK(isamek(qirr, qbz(:,iq_bz), g0), "isamek")
230 122 : qirredtofull(iq_ibz) = iq_bz
231 : else
232 0 : ABI_ERROR(sjoin("Full BZ does not contain IBZ q-point:", ktoa(qirr)))
233 : end if
234 : end do
235 :
236 : ! Build qpttoqpt table. See also mkqptequiv
237 40 : ABI_MALLOC(qpttoqpt, (2, cryst%nsym, nqbz))
238 222730 : qpttoqpt = -1
239 1546 : do iq_bz=1,nqbz
240 75274 : do isym=1,cryst%nsym
241 1843200 : tmp_qpt = matmul(cryst%symrec(:,:,isym), qbz(:,iq_bz))
242 :
243 73728 : isq_bz = qrank%get_index(tmp_qpt)
244 73728 : if (isq_bz == -1) then
245 0 : ABI_ERROR("Looks like no kpoint equiv to q by symmetry without time reversal!")
246 : end if
247 73728 : qpttoqpt(1,isym,isq_bz) = iq_bz
248 :
249 : ! q --> -q
250 294912 : tmp_qpt = -tmp_qpt
251 73728 : isq_bz = qrank%get_index(tmp_qpt)
252 73728 : if (isq_bz == -1) then
253 0 : ABI_ERROR("Looks like no kpoint equiv to q by symmetry with time reversal!")
254 : end if
255 75264 : qpttoqpt(2,isym,isq_bz) = iq_bz
256 : end do
257 : end do
258 :
259 10 : call qrank%free()
260 :
261 10 : end subroutine ephtk_mkqtabs
262 : !!***
263 :
264 : !----------------------------------------------------------------------
265 :
266 : !!****f* m_ephtk/ephtk_gam_atm2qnu
267 : !! NAME
268 : !! ephtk_gam_atm2qnu
269 : !!
270 : !! FUNCTION
271 : !! This routine takes the gamma matrices in the atomic representation and
272 : !! multiplies them by the displ_red matrices. Based on gam_mult_displ
273 : !!
274 : !! INPUTS
275 : !! natom3 = number of phonon branches (3*natom)
276 : !! displ_red = phonon mode displacement vectors in reduced coordinates.
277 : !! gam_bare = bare gamma matrices before multiplication
278 : !!
279 : !! OUTPUT
280 : !! gam_now = output gamma matrices multiplied by displacement matrices
281 : !!
282 : !! SOURCE
283 :
284 15200 : subroutine ephtk_gam_atm2qnu(natom3, displ_red, gam_atm, gam_qnu)
285 :
286 : !Arguments -------------------------------
287 : integer, intent(in) :: natom3
288 : real(dp), intent(in) :: displ_red(2,natom3,natom3), gam_atm(2,natom3,natom3)
289 : real(dp), intent(out) :: gam_qnu(natom3)
290 :
291 : !Local variables -------------------------
292 : integer,save :: enough = 0
293 : integer :: nu
294 : character(len=500) :: msg
295 30400 : real(dp) :: zgemm_tmp_mat(2,natom3,natom3), gam_now(2,natom3,natom3)
296 : ! *********************************************************************
297 :
298 15200 : call zgemm('c','n',natom3, natom3, natom3, cone, displ_red, natom3, gam_atm, natom3, czero, zgemm_tmp_mat, natom3)
299 :
300 471200 : gam_now = zero
301 15200 : call zgemm('n','n',natom3,natom3,natom3,cone,zgemm_tmp_mat,natom3,displ_red,natom3,czero,gam_now,natom3)
302 :
303 : ! Extract gamma(q,nu)
304 60800 : do nu=1,natom3
305 45600 : gam_qnu(nu) = gam_now(1, nu, nu)
306 60800 : if (abs(gam_now(2, nu, nu)) > tol8) then
307 0 : enough = enough + 1
308 0 : if (enough <= 30) then
309 0 : write (msg,'(a,i0,a,es16.8)')' non-zero imaginary part for branch: ',nu,', img: ',gam_now(2, nu, nu)
310 0 : ABI_WARNING(msg)
311 : end if
312 : end if
313 : end do
314 :
315 15200 : end subroutine ephtk_gam_atm2qnu
316 : !!***
317 :
318 : !----------------------------------------------------------------------
319 :
320 : !!****f* m_ephtk/ephtk_gkknu_from_atm
321 : !! NAME
322 : !! ephtk_gkknu_from_atm
323 : !!
324 : !! FUNCTION
325 : !! Transform the gkk matrix elements from (atom, red_direction) basis to phonon-mode basis.
326 : !!
327 : !! INPUTS
328 : !! nb1,nb2=Number of bands in gkq_atm matrix.
329 : !! nk=Number of k-points (usually 1)
330 : !! natom=Number of atoms.
331 : !! gkq_atm(2,nb1,nb2,3*natom)=EPH matrix elements in the atomic basis.
332 : !! phfrq(3*natom)=Phonon frequencies in Ha
333 : !! displ_red(2,3*natom,3*natom)=Phonon displacement in reduced coordinates.
334 : !!
335 : !! OUTPUT
336 : !! gkq_nu(2,nb1,nb2,3*natom)=EPH matrix elements in the phonon-mode basis.
337 : !!
338 : !! SOURCE
339 :
340 135493 : subroutine ephtk_gkknu_from_atm(nb1, nb2, nk, natom, gkq_atm, phfrq, displ_red, gkq_nu)
341 :
342 : !Arguments ------------------------------------
343 : !scalars
344 : integer,intent(in) :: nb1, nb2, nk, natom
345 : !arrays
346 : real(dp),intent(in) :: phfrq(3*natom), displ_red(2,3*natom,3*natom)
347 : real(dp),intent(in) :: gkq_atm(2,nb1,nb2,nk,3*natom)
348 : real(dp),intent(out) :: gkq_nu(2,nb1,nb2,nk,3*natom)
349 :
350 : !Local variables-------------------------
351 : integer :: nu,ipc
352 : ! *************************************************************************
353 :
354 16285561 : gkq_nu = zero
355 :
356 : ! Loop over phonon branches.
357 945895 : do nu=1,3*natom
358 : ! Ignore negative or too small frequencies
359 810402 : if (phfrq(nu) < EPHTK_WTOL) cycle
360 :
361 : ! Transform the gkk from (atom, reduced direction) basis to phonon mode representation.
362 5599368 : do ipc=1,3*natom
363 : gkq_nu(1,:,:,:,nu) = gkq_nu(1,:,:,:,nu) &
364 : + gkq_atm(1,:,:,:,ipc) * displ_red(1,ipc,nu) &
365 46034268 : - gkq_atm(2,:,:,:,ipc) * displ_red(2,ipc,nu)
366 : gkq_nu(2,:,:,:,nu) = gkq_nu(2,:,:,:,nu) &
367 : + gkq_atm(1,:,:,:,ipc) * displ_red(2,ipc,nu) &
368 46834722 : + gkq_atm(2,:,:,:,ipc) * displ_red(1,ipc,nu)
369 : end do
370 :
371 16062697 : gkq_nu(:,:,:,:,nu) = gkq_nu(:,:,:,:,nu) / sqrt(two * phfrq(nu))
372 : end do
373 :
374 135493 : end subroutine ephtk_gkknu_from_atm
375 : !!***
376 :
377 : !----------------------------------------------------------------------
378 :
379 : !!****f* m_ephtk/ephtk_update_ebands
380 : !! NAME
381 : !! ephtk_update_ebands
382 : !!
383 : !! FUNCTION
384 : !! Update ebands according to dtset%occopt, tsmear, mbpt_sciss, eph_fermie, eph_extrael
385 : !!
386 : !! INPUTS
387 : !! dtset<dataset_type>=All input variables for this dataset.
388 : !! fileqpdatain: QPDATA file
389 : !!
390 : !! SOURCE
391 :
392 99 : subroutine ephtk_update_ebands(dtset, ebands, header)
393 :
394 : !Arguments ------------------------------------
395 : !scalars
396 : type(dataset_type),intent(in) :: dtset
397 : type(ebands_t),intent(inout) :: ebands
398 : character(len=*),intent(in) :: header
399 :
400 : !Local variables-------------------------
401 : !scalars
402 : real(dp),parameter :: nholes = zero
403 : character(len=500) :: msg
404 : integer :: units(2)
405 : ! *************************************************************************
406 :
407 297 : units = [std_out, ab_out]
408 :
409 99 : if (abs(dtset%mbpt_sciss) > tol6) then
410 : ! Apply the scissor operator
411 : call wrtout(units, sjoin(" Applying scissors operator to the conduction states with value: ", &
412 0 : ftoa(dtset%mbpt_sciss * Ha_eV, fmt="(f6.2)"), " (eV)"))
413 0 : call ebands%apply_scissors(dtset%mbpt_sciss)
414 : end if
415 :
416 99 : if (dtset%occopt /= ebands%occopt .or. abs(dtset%tsmear - ebands%tsmear) > tol12) then
417 : write(msg,"(2a,2(a,i0,a,f14.6,a))")&
418 21 : " Changing occupation scheme as input occopt and tsmear differ from those read from WFK file.",ch10,&
419 21 : " From WFK file: occopt = ",ebands%occopt,", tsmear = ",ebands%tsmear,ch10,&
420 42 : " From input: occopt = ",dtset%occopt,", tsmear = ",dtset%tsmear,ch10
421 21 : call wrtout(units, msg)
422 21 : call ebands%set_scheme(dtset%occopt, dtset%tsmear, dtset%spinmagntarget, dtset%prtvol)
423 : end if
424 :
425 : ! Default value of eph_fermie is zero hence no tolerance is used!
426 99 : if (dtset%eph_fermie /= zero) then
427 1 : ABI_CHECK(dtset%eph_extrael == zero, "eph_fermie and eph_extrael are mutually exclusive")
428 1 : call wrtout(units, sjoin(" Fermi level set by the user at:", ftoa(dtset%eph_fermie * Ha_eV, fmt="(f6.2)"), " (eV)"))
429 1 : call ebands%set_fermie(dtset%eph_fermie, msg)
430 1 : call wrtout(units, msg)
431 :
432 98 : else if (abs(dtset%eph_extrael) > zero) then
433 20 : call wrtout(units, sjoin(" Adding eph_extrael:", ftoa(dtset%eph_extrael), "to input nelect:", ftoa(ebands%nelect)))
434 20 : call ebands%set_scheme(dtset%occopt, dtset%tsmear, dtset%spinmagntarget, dtset%prtvol, update_occ=.False.)
435 20 : call ebands%set_extrael(dtset%eph_extrael, nholes, dtset%spinmagntarget, msg)
436 20 : call wrtout(units, msg)
437 : end if
438 :
439 : ! Recompute occupations. This is needed if WFK files have been produced in a NSCF run
440 : ! since occ are set to zero, and fermie is taken from the previous density.
441 99 : if (dtset%kptopt > 0) then
442 98 : call ebands%update_occ(dtset%spinmagntarget, prtvol=dtset%prtvol)
443 196 : call ebands%print([std_out], header=header, prtvol=dtset%prtvol)
444 : end if
445 :
446 99 : end subroutine ephtk_update_ebands
447 : !!***
448 :
449 : !----------------------------------------------------------------------
450 :
451 : !!****f* m_ephtk/ephtk_get_mpw_gmax
452 : !! NAME
453 : !! ephtk_get_mpw_gmax
454 : !!
455 : !! FUNCTION
456 : !! Compute maximum number of plane-waves over k and k+q where k and k+q are in the BZ.
457 : !! we also need the max components of the G-spheres (k, k+q) in order to allocate the workspace array work
458 : !! used to symmetrize the wavefunctions in G-space.
459 : !! Note that we loop over the full BZ instead of the IBZ(k)
460 : !! This part is slow for very dense meshes, should try to use a geometrical approach...
461 : !!
462 : !! INPUTS
463 : !!
464 : !! SOURCE
465 :
466 60 : subroutine ephtk_get_mpw_gmax(nkpt, kpts, ecut, gmet, mpw, gmax, comm, init_with_zero)
467 :
468 : !Arguments ------------------------------------
469 : integer,intent(in) :: nkpt
470 : integer,intent(out) :: mpw, gmax(3)
471 : real(dp),intent(in) :: ecut, kpts(3,nkpt), gmet(3,3)
472 : integer,intent(in) :: comm
473 : logical,optional,intent(in) :: init_with_zero
474 :
475 : !Local variables ------------------------------
476 : integer,parameter :: istwfk1 = 1
477 : integer :: ik,i1,i2,i3,cnt,ipw,ii,onpw,my_mpw,my_gmax(3),ierr, my_rank, nprocs
478 : real(dp) :: kk(3), kq(3)
479 60 : integer,allocatable :: gtmp(:,:)
480 : logical :: init_with_zero__
481 : ! *************************************************************************
482 :
483 60 : my_rank = xmpi_comm_rank(comm); nprocs = xmpi_comm_size(comm)
484 :
485 60 : init_with_zero__ = .True.; if (present(init_with_zero)) init_with_zero__ = init_with_zero
486 :
487 2 : if (init_with_zero__) then
488 60 : mpw = 0; gmax = 0
489 : end if
490 :
491 : cnt = 0
492 358 : do ik=1,nkpt
493 1192 : kk = kpts(:, ik)
494 1252 : do i3=-1,1
495 3874 : do i2=-1,1
496 11622 : do i1=-1,1
497 8046 : cnt = cnt + 1; if (mod(cnt, nprocs) /= my_rank) cycle ! MPI parallelism inside comm
498 32184 : kq = kk + half * [i1, i2, i3]
499 : ! TODO: g0 umklapp here can enter into play gmax may not be large enough!
500 8046 : call get_kg(kq, istwfk1, 1.1_dp * ecut, gmet, onpw, gtmp)
501 8046 : mpw = max(mpw, onpw)
502 6937057 : do ipw=1,onpw
503 27724090 : do ii=1,3
504 27716044 : gmax(ii) = max(gmax(ii), abs(gtmp(ii, ipw)))
505 : end do
506 : end do
507 10728 : ABI_FREE(gtmp)
508 : end do
509 : end do
510 : end do
511 : end do
512 :
513 60 : my_mpw = mpw; call xmpi_max(my_mpw, mpw, comm, ierr)
514 60 : my_gmax = gmax; call xmpi_max(my_gmax, gmax, comm, ierr)
515 :
516 60 : end subroutine ephtk_get_mpw_gmax
517 : !!***
518 :
519 : !!****f* m_epthk/ephtk_v1atm_to_vqnu
520 : !! NAME
521 : !! ephtk_v1atm_to_vqnu
522 : !!
523 : !! FUNCTION
524 : !! Receive potentials in atomic representation and return potential in phonon representation
525 : !!
526 : !! INPUTS
527 : !!
528 : !! OUTPUT
529 : !!
530 : !! SOURCE
531 :
532 0 : pure subroutine ephtk_v1atm_to_vqnu(cplex, nfft, nspden, natom3, v1_atm, displ_red, v1_qnu)
533 :
534 : !Arguments ------------------------------------
535 : !scalars
536 : integer,intent(in) :: cplex, nfft, nspden, natom3
537 : !arrays
538 : real(dp),intent(in) :: v1_atm(cplex, nfft, nspden, natom3)
539 : real(dp),intent(out) :: v1_qnu(2, nfft, nspden, natom3)
540 : real(dp),intent(in) :: displ_red(2, natom3, natom3)
541 :
542 : !Local variables-------------------------------
543 : !scalars
544 : integer :: nu, ip, ispden
545 : !************************************************************************
546 :
547 0 : do nu=1,natom3
548 : ! v1_qnu = \sum_{ka} phdispl{ka}(q,nu) D_{ka,q} V_scf(r)
549 : ! NOTE: prefactor 1/sqrt(2 w(q,nu)) is not included in the potentials.
550 : ! v1_qnu(2, nfft, nspden, natom3), v1_atm(cplex, nfft, nspden, natom3)
551 0 : v1_qnu(:, :, :, nu) = zero
552 0 : do ip=1,natom3
553 0 : do ispden=1,nspden
554 0 : if (cplex == 2) then
555 : v1_qnu(1, :, ispden, nu) = v1_qnu(1, :, ispden, nu) + &
556 0 : displ_red(1,ip,nu) * v1_atm(1,:,ispden,ip) - displ_red(2,ip,nu) * v1_atm(2,:,ispden,ip)
557 : v1_qnu(2, :, ispden, nu) = v1_qnu(2, :, ispden, nu) + &
558 0 : displ_red(2,ip,nu) * v1_atm(1,:,ispden,ip) + displ_red(1,ip,nu) * v1_atm(2,:,ispden,ip)
559 : else
560 : ! Gamma point. d(q) = d(-q)* --> d is real.
561 0 : v1_qnu(1, :, ispden, nu) = v1_qnu(1, :, ispden, nu) + displ_red(1,ip,nu) * v1_atm(1,:,ispden,ip)
562 : end if
563 : end do
564 : end do
565 : end do
566 :
567 0 : end subroutine ephtk_v1atm_to_vqnu
568 : !!***
569 :
570 : !!****f* m_ephtk/ephtk_skip_phmode
571 : !! NAME
572 : !! ephtk_skip_mode
573 : !!
574 : !! FUNCTION
575 : !! Ignore contribution of phonon mode depending on phonon frequency value or mode index.
576 : !!
577 : !! INPUTS
578 : !! nu: mode index
579 : !! wqnu: phonon frequency
580 : !! eph_phrange_w: range for phonon frequency.
581 : !!
582 : !! SOURCE
583 :
584 155448 : pure logical function ephtk_skip_phmode(nu, wqnu, phmodes_skip, eph_phrange_w) result(skip)
585 :
586 : !Arguments ------------------------------------
587 : integer,intent(in) :: nu, phmodes_skip(:)
588 : real(dp),intent(in) :: wqnu, eph_phrange_w(2)
589 : ! *************************************************************************
590 :
591 155448 : skip = wqnu < EPHTK_WTOL .or. phmodes_skip(nu) == 1
592 :
593 : ! Check frequency range
594 155448 : if (abs(eph_phrange_w(2)) > tol12) then
595 1716 : if (eph_phrange_w(2) > zero) then
596 : ! wqnu must be inside range
597 1716 : skip = skip .or. .not. (wqnu >= eph_phrange_w(1) .and. wqnu <= eph_phrange_w(2))
598 : else
599 : ! wqnu must be outside range
600 0 : skip = skip .or. (wqnu >= eph_phrange_w(1) .and. wqnu <= eph_phrange_w(2))
601 : end if
602 : end if
603 :
604 155448 : end function ephtk_skip_phmode
605 : !!***
606 :
607 73728 : end module m_ephtk
608 : !!***
|
