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