Line data Source code
1 : !!****m* ABINIT/m_gkk
2 : !! NAME
3 : !!
4 : !! FUNCTION
5 : !! Tools for the computation of electron-phonon coupling matrix elements (gkk)
6 : !!
7 : !! COPYRIGHT
8 : !! Copyright (C) 2008-2026 ABINIT group (GKA, MG)
9 : !! This file is distributed under the terms of the
10 : !! GNU General Public License, see ~abinit/COPYING
11 : !! or http://www.gnu.org/copyleft/gpl.txt .
12 : !!
13 : !! SOURCE
14 :
15 : #if defined HAVE_CONFIG_H
16 : #include "config.h"
17 : #endif
18 :
19 : #include "abi_common.h"
20 :
21 : module m_gkk
22 :
23 : use defs_basis
24 : use m_abicore
25 : use m_xmpi
26 : use m_errors
27 : use m_dtset
28 : use m_ifc
29 : use m_ebands
30 : use m_ddb
31 : use m_dvdb
32 : use m_fft
33 : use m_hamiltonian
34 : use m_pawcprj
35 : use m_wfk
36 : use m_nctk
37 : use m_dtfil
38 : use netcdf
39 :
40 : use defs_abitypes, only : MPI_type
41 : use m_time, only : cwtime, sec2str
42 : use m_io_tools, only : iomode_from_fname
43 : use m_fstrings, only : itoa, sjoin, ktoa, ltoa, strcat
44 : use m_symtk, only : littlegroup_q
45 : use m_fftcore, only : get_kg
46 : use defs_datatypes, only : pseudopotential_type
47 : use m_crystal, only : crystal_t
48 : use m_bz_mesh, only : findqg0
49 : use m_cgtools, only : dotprod_g
50 : use m_kg, only : getph
51 : use m_pawang, only : pawang_type
52 : use m_pawrad, only : pawrad_type
53 : use m_pawtab, only : pawtab_type
54 : use m_pawfgr, only : pawfgr_type
55 : use m_eig2d, only : gkk_t, gkk_init, gkk_ncwrite, gkk_free
56 : use m_wfd, only : wfd_t
57 : use m_getgh1c, only : getgh1c, rf_transgrid_and_pack, getgh1c_setup
58 : use m_ephtk, only : ephtk_v1atm_to_vqnu
59 :
60 : implicit none
61 :
62 : private
63 : !!***
64 :
65 : public :: eph_gkk
66 : public :: ncwrite_v1qnu ! Compute \delta V_{q,nu)(r) and dump results to netcdf file.
67 :
68 : contains !===========================================================================
69 : !!***
70 :
71 : !!****f* m_gkk/eph_gkk
72 : !! NAME
73 : !! eph_gkk
74 : !!
75 : !! FUNCTION
76 : !! Compute electron-phonon coupling matrix elements.
77 : !!
78 : !! INPUTS
79 : !! wk0_path=String with the path to the GS unperturbed WFK file.
80 : !! ngfft(18),ngfftf(18)=Coarse and Fine FFT meshes.
81 : !! dtset<dataset_type>=All input variables for this dataset.
82 : !! ebands<ebands_t>=The GS KS band structure (energies, occupancies, k-weights...)
83 : !! dvdb<dbdb_type>=Database with the DFPT SCF potentials.
84 : !! ifc<ifc_type>=interatomic force constants and corresponding real space grid info.
85 : !! pawfgr <type(pawfgr_type)>=fine grid parameters and related data
86 : !! pawang<pawang_type)>=PAW angular mesh and related data.
87 : !! pawrad(ntypat*usepaw)<pawrad_type>=Paw radial mesh and related data.
88 : !! pawtab(ntypat*usepaw)<pawtab_type>=Paw tabulated starting data.
89 : !! psps<pseudopotential_type>=Variables related to pseudopotentials.
90 : !! comm=MPI communicator.
91 : !!
92 : !! OUTPUT
93 : !!
94 : !! SOURCE
95 :
96 2 : subroutine eph_gkk(wfk0_path,wfq_path,dtfil,ngfft,ngfftf,dtset,cryst,ebands_k,ebands_kq,dvdb,ifc,&
97 2 : pawfgr,pawang,pawrad,pawtab,psps,mpi_enreg,comm)
98 :
99 : !Arguments ------------------------------------
100 : !scalars
101 : character(len=*),intent(in) :: wfk0_path, wfq_path
102 : integer,intent(in) :: comm
103 : type(datafiles_type),intent(in) :: dtfil
104 : type(dataset_type),intent(in) :: dtset
105 : type(crystal_t),intent(in) :: cryst
106 : type(ebands_t),intent(in) :: ebands_k, ebands_kq
107 : type(dvdb_t),target,intent(inout) :: dvdb
108 : type(pawang_type),intent(in) :: pawang
109 : type(pseudopotential_type),intent(in) :: psps
110 : type(pawfgr_type),intent(in) :: pawfgr
111 : type(ifc_type),intent(in) :: ifc
112 : type(mpi_type),intent(inout) :: mpi_enreg
113 : !arrays
114 : integer,intent(in) :: ngfft(18),ngfftf(18)
115 : type(pawrad_type),intent(in) :: pawrad(psps%ntypat*psps%usepaw)
116 : type(pawtab_type),intent(in) :: pawtab(psps%ntypat*psps%usepaw)
117 :
118 : !Local variables ------------------------------
119 : !scalars
120 : integer,parameter :: tim_getgh1c=1, berryopt0=0, useylmgr1=0, master=0, qptopt1 = 1
121 : integer :: my_rank,nproc,mband,mband_kq,my_minb,my_maxb,nsppol,nkpt,nkpt_kq,idir,ipert
122 : integer :: cplex,db_iqpt,natom,natom3,ipc,nspinor
123 : integer :: ib1,ib2,band,ik,ikq,timerev_q
124 : integer :: spin,istwf_k,istwf_kq,npw_k,npw_kq, comm_rpt
125 : integer :: mpw,mpw_k,mpw_kq,ierr,my_kstart,my_kstop,ncid
126 : integer :: n1,n2,n3,n4,n5,n6,nspden,ncerr
127 : integer :: sij_opt,usecprj,usevnl,optlocal,optnl,opt_gvnlx1
128 : integer :: nfft,nfftf,mgfft,mgfftf,nkpg,nkpg1, interpolated
129 : real(dp) :: cpu,wall,gflops,ecut,eshift,eig0nk,dotr,doti
130 : logical :: i_am_master, gen_eigenpb
131 2 : type(wfd_t) :: wfd_k, wfd_kq
132 2 : type(gs_hamiltonian_type) :: gs_hamkq
133 2 : type(rf_hamiltonian_type) :: rf_hamkq
134 2 : type(gkk_t) :: gkk2d
135 : character(len=500) :: msg, what
136 : character(len=fnlen) :: fname, gkkfilnam
137 : !arrays
138 4 : integer :: g0_k(3),symq(4,2,cryst%nsym), units(2)
139 2 : integer,allocatable :: kg_k(:,:),kg_kq(:,:),nband(:,:),nband_kq(:,:),wfd_istwfk(:)
140 : real(dp) :: ylmgr_kq_dum(1,1,1) ! ylmgr_k_dum(1,1,1),
141 4 : real(dp) :: kk(3),kq(3),qpt(3),phfrq(3*cryst%natom),dvdb_qdamp(1)
142 2 : real(dp),allocatable :: displ_cart(:,:,:),displ_red(:,:,:), eigens_kq(:,:,:)
143 2 : real(dp),allocatable :: grad_berry(:,:),kinpw_kq(:),kpg_kq(:,:),kpg_k(:,:),dkinpw(:)
144 2 : real(dp),allocatable :: ffnl_k(:,:,:,:),ffnl_kq(:,:,:,:),ph3d_k(:,:,:),ph3d_kq(:,:,:)
145 2 : real(dp),allocatable :: v1scf(:,:,:,:),gkk(:,:,:,:,:), bras(:,:,:),kets(:,:,:),h1_kets(:,:,:)
146 2 : real(dp),allocatable :: ph1d(:,:),vlocal(:,:,:,:),vlocal1(:,:,:,:,:), ylm_kq(:,:),ylm_k(:,:)
147 2 : real(dp),allocatable :: dummy_vtrial(:,:),gvnlx1(:,:), gs1c(:,:), gkq_atm(:,:,:,:)
148 2 : logical,allocatable :: bks_mask(:,:,:),bks_mask_kq(:,:,:),keep_ur(:,:,:),keep_ur_kq(:,:,:)
149 2 : type(pawcprj_type),allocatable :: cwaveprj0(:,:) !natom,nspinor*usecprj)
150 : !************************************************************************
151 :
152 6 : units = [std_out, ab_out]
153 :
154 2 : what = "(GKK files)"; if (dtset%eph_task == -2) what = "GKQ file"
155 2 : write(msg, '(3a)') " Computation of electron-phonon coupling matrix elements ", trim(what), ch10
156 2 : call wrtout(units, msg, do_flush=.True.)
157 :
158 2 : if (psps%usepaw == 1) then
159 0 : ABI_ERROR("PAW not implemented")
160 : ABI_UNUSED((/pawang%nsym, pawrad(1)%mesh_size/))
161 : end if
162 :
163 2 : my_rank = xmpi_comm_rank(comm); nproc = xmpi_comm_size(comm); i_am_master = my_rank == master
164 :
165 : ! Copy important dimensions
166 2 : natom = cryst%natom; natom3 = 3 * natom; nsppol = ebands_k%nsppol; nspinor = ebands_k%nspinor; nspden = dtset%nspden
167 2 : nkpt = ebands_k%nkpt; mband = ebands_k%mband; nkpt_kq = ebands_kq%nkpt; mband_kq = ebands_kq%mband
168 2 : ecut = dtset%ecut
169 : !write(std_out, *)"ebands dims (b, k, s): ", ebands_k%mband, ebands_k%nkpt, ebands_k%nsppol
170 : !write(std_out, *)"ebands_kq dims (b, k, s): ", ebands_kq%mband, ebands_kq%nkpt, ebands_kq%nsppol
171 :
172 8 : qpt = dtset%qptn(:)
173 :
174 8 : nfftf = product(ngfftf(1:3)); mgfftf = maxval(ngfftf(1:3))
175 14 : nfft = product(ngfft(1:3)) ; mgfft = maxval(ngfft(1:3))
176 2 : n1=ngfft(1); n2=ngfft(2); n3=ngfft(3)
177 2 : n4=ngfft(4); n5=ngfft(5); n6=ngfft(6)
178 :
179 : ! Open the DVDB file
180 2 : call dvdb%open_read(ngfftf, xmpi_comm_self)
181 :
182 : ! Initialize the wave function descriptors.
183 : ! For the time being, no memory distribution, each node has the full set of states.
184 2 : my_minb = 1; my_maxb = mband
185 :
186 8 : ABI_MALLOC(nband, (nkpt, nsppol))
187 10 : ABI_MALLOC(bks_mask,(mband, nkpt, nsppol))
188 8 : ABI_MALLOC(keep_ur,(mband, nkpt ,nsppol))
189 216 : nband=mband; bks_mask=.False.; keep_ur=.False.
190 :
191 8 : ABI_MALLOC(nband_kq, (nkpt_kq, nsppol))
192 10 : ABI_MALLOC(bks_mask_kq,(mband_kq, nkpt_kq, nsppol))
193 8 : ABI_MALLOC(keep_ur_kq,(mband_kq, nkpt_kq ,nsppol))
194 216 : nband_kq=mband_kq; bks_mask_kq=.False.; keep_ur_kq=.False.
195 :
196 : ! Distribute the k-points over the processors
197 2 : call xmpi_split_work(nkpt,comm,my_kstart,my_kstop)
198 18 : do ik=1,nkpt
199 16 : if (.not. (ik >= my_kstart .and. ik <= my_kstop)) cycle
200 64 : kk = ebands_k%kptns(:,ik)
201 64 : kq = kk + qpt
202 : ! Find the index of the k+q point
203 16 : call findqg0(ikq,g0_k,kq,nkpt_kq,ebands_kq%kptns(:,:), [1,1,1])
204 112 : bks_mask(:,ik,:) = .True.
205 114 : bks_mask_kq(:,ikq,:) = .True.
206 : end do
207 :
208 : ! Impose istwfk=1 for all k points. This is also done in respfn (see inkpts)
209 : ! wfd_read_wfk will handle a possible conversion if WFK contains istwfk /= 1.
210 6 : ABI_MALLOC(wfd_istwfk, (nkpt))
211 18 : wfd_istwfk = 1
212 :
213 : ! Initialize the wavefunction descriptors
214 : call wfd_k%init(cryst,pawtab,psps,keep_ur,mband,nband,nkpt,nsppol,bks_mask,&
215 : nspden,nspinor,ecut,dtset%ecutsm,dtset%dilatmx,wfd_istwfk,ebands_k%kptns,ngfft,&
216 2 : dtset%nloalg,dtset%prtvol,dtset%pawprtvol,comm)
217 2 : ABI_FREE(wfd_istwfk)
218 :
219 4 : call wfd_k%print([std_out], header="Wavefunctions on the k-points grid")
220 :
221 6 : ABI_MALLOC(wfd_istwfk, (nkpt_kq))
222 18 : wfd_istwfk = 1
223 :
224 : call wfd_kq%init(cryst,pawtab,psps,keep_ur_kq,mband_kq,nband_kq,nkpt_kq,nsppol,bks_mask_kq,&
225 : nspden,nspinor,ecut,dtset%ecutsm,dtset%dilatmx,wfd_istwfk,ebands_kq%kptns,ngfft,&
226 2 : dtset%nloalg,dtset%prtvol,dtset%pawprtvol,comm)
227 2 : ABI_FREE(wfd_istwfk)
228 :
229 4 : call wfd_kq%print([std_out], header="Wavefunctions on the q-shifted k-points grid")
230 :
231 2 : ABI_FREE(nband)
232 2 : ABI_FREE(bks_mask)
233 2 : ABI_FREE(keep_ur)
234 2 : ABI_FREE(nband_kq)
235 2 : ABI_FREE(bks_mask_kq)
236 2 : ABI_FREE(keep_ur_kq)
237 :
238 : ! Read wavefunctions on the k-points grid and q-shifted k-points grid.
239 2 : call wfd_k%read_wfk(wfk0_path, iomode_from_fname(wfk0_path))
240 2 : call wfd_kq%read_wfk(wfq_path, iomode_from_fname(wfq_path))
241 :
242 : ! ph1d(2,3*(2*mgfft+1)*natom)=one-dimensional structure factor information on the coarse grid.
243 6 : ABI_MALLOC(ph1d, (2,3*(2*mgfft+1)*natom))
244 2 : call getph(cryst%atindx,natom,n1,n2,n3,ph1d,cryst%xred)
245 :
246 : ! Find the appropriate value of mpw
247 2 : call find_mpw(mpw_k, ebands_k%kptns(:,:), nsppol, nkpt, cryst%gmet,ecut,comm)
248 2 : call find_mpw(mpw_kq, ebands_kq%kptns(:,:), nsppol, nkpt_kq, cryst%gmet,ecut,comm)
249 2 : mpw = max(mpw_k, mpw_kq)
250 :
251 : ! Allow PW-arrays dimensioned with mpw
252 6 : ABI_MALLOC(kg_k, (3, mpw))
253 4 : ABI_MALLOC(kg_kq, (3, mpw))
254 :
255 : ! Spherical Harmonics for useylm==1.
256 : ! TODO: These arrays should be allocated with npw_k and npw_kq
257 8 : ABI_MALLOC(ylm_k,(mpw, psps%mpsang*psps%mpsang*psps%useylm))
258 6 : ABI_MALLOC(ylm_kq,(mpw, psps%mpsang*psps%mpsang*psps%useylm))
259 :
260 : ! TODO FOR PAW
261 2 : usecprj = 0
262 4 : ABI_MALLOC(cwaveprj0, (natom, nspinor*usecprj))
263 :
264 : ! Prepare call to getgh1c
265 2 : usevnl = 0
266 2 : optlocal = 1 ! local part of H^(1) is computed in gh1c=<G|H^(1)|C>
267 2 : optnl = 2 ! non-local part of H^(1) is totally computed in gh1c=<G|H^(1)|C>
268 2 : opt_gvnlx1 = 0 ! gvnlx1 is output
269 2 : ABI_MALLOC(gvnlx1, (2,usevnl))
270 2 : ABI_MALLOC(grad_berry, (2,nspinor*(berryopt0/4)))
271 :
272 : ! This part is taken from dfpt_vtorho
273 : !==== Initialize most of the Hamiltonian (and derivative) ====
274 : !1) Allocate all arrays and initialize quantities that do not depend on k and spin.
275 : !2) Perform the setup needed for the non-local factors:
276 : !* Norm-conserving: Constant kleimann-Bylander energies are copied from psps to gs_hamk.
277 : !* PAW: Initialize the overlap coefficients and allocate the Dij coefficients.
278 :
279 : call gs_hamkq%init(psps,pawtab,nspinor,nsppol,nspden,natom,&
280 : dtset%typat,cryst%xred,nfft,mgfft,ngfft,cryst%rprimd,dtset%nloalg,&
281 : usecprj=usecprj,ph1d=ph1d,nucdipmom=dtset%nucdipmom,gpu_option=dtset%gpu_option,&
282 2 : comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab,mpi_spintab=mpi_enreg%my_isppoltab)
283 :
284 : ! Allocate vlocal. Note nvloc
285 12 : ABI_MALLOC(vlocal,(n4,n5,n6,gs_hamkq%nvloc))
286 : ! Allocate work space arrays.
287 3466 : ABI_CALLOC(dummy_vtrial, (nfftf,nspden))
288 :
289 2 : call cwtime(cpu, wall, gflops, "start")
290 :
291 2 : interpolated = 0
292 2 : if (dtset%eph_use_ftinterp /= 0) then
293 0 : ABI_WARNING(sjoin("Enforcing FT interpolation for q-point", ktoa(qpt)))
294 0 : comm_rpt = xmpi_comm_self
295 0 : call dvdb%ftinterp_setup(dtset%ddb_ngqpt, qptopt1, 1, dtset%ddb_shiftq, nfftf, ngfftf, comm_rpt)
296 0 : cplex = 2
297 0 : ABI_MALLOC(v1scf, (cplex, nfftf, nspden, dvdb%my_npert))
298 0 : call dvdb%ftinterp_qpt(qpt, nfftf, ngfftf, v1scf, dvdb%comm_rpt)
299 0 : interpolated = 1
300 : else
301 : ! Find the index of the q-point in the DVDB.
302 2 : db_iqpt = dvdb%findq(qpt)
303 2 : if (db_iqpt /= -1) then
304 2 : if (dtset%prtvol > 0) call wrtout(std_out, sjoin("Found: ",ktoa(qpt)," in DVDB with index ",itoa(db_iqpt)))
305 : ! Read or reconstruct the dvscf potentials for all 3*natom perturbations.
306 : ! This call allocates v1scf(cplex, nfftf, nspden, 3*natom))
307 2 : call dvdb%readsym_allv1(db_iqpt, cplex, nfftf, ngfftf, v1scf, comm)
308 : else
309 0 : ABI_WARNING(sjoin("Cannot find q-point:", ktoa(qpt), "in DVDB file"))
310 : end if
311 : end if
312 :
313 : ! Examine the symmetries of the q wavevector
314 2 : call littlegroup_q(cryst%nsym,qpt,symq,cryst%symrec,cryst%symafm,timerev_q,prtvol=dtset%prtvol)
315 :
316 : ! Allocate vlocal1 with correct cplex. Note nvloc
317 14 : ABI_MALLOC_OR_DIE(vlocal1,(cplex*n4,n5,n6,gs_hamkq%nvloc,natom3), ierr)
318 :
319 10 : ABI_MALLOC(displ_cart, (2,3*cryst%natom,3*cryst%natom))
320 6 : ABI_MALLOC(displ_red, (2,3*cryst%natom,3*cryst%natom))
321 :
322 2 : if (dtset%eph_task == 2) then
323 : ! Write GKK files (1 file for perturbation)
324 5 : ABI_MALLOC(gkk, (2*mband*nsppol,nkpt,1,1,mband_kq))
325 :
326 1 : else if (dtset%eph_task == -2) then
327 : ! Write GKQ file with all perturbations. gkq are given in the atom representation.
328 : ! TODO: Assuming mband_kq == mband
329 5 : ABI_MALLOC(gkq_atm, (2, mband_kq, mband, nkpt))
330 1 : if (i_am_master) then
331 1 : call ifc%fourq(cryst, qpt, phfrq, displ_cart, out_displ_red=displ_red)
332 1 : fname = strcat(dtfil%filnam_ds(4), "_GKQ.nc")
333 1 : NCF_CHECK_MSG(nctk_open_create(ncid, fname, xmpi_comm_self), "Creating GKQ file")
334 1 : NCF_CHECK(cryst%ncwrite(ncid))
335 : ! Write bands on k mesh.
336 1 : NCF_CHECK(ebands_k%ncwrite(ncid))
337 2 : ncerr = nctk_def_dims(ncid, [nctkdim_t('number_of_phonon_modes', natom3)], defmode=.True.)
338 1 : NCF_CHECK(ncerr)
339 : ncerr = nctk_def_iscalars(ncid, [character(len=nctk_slen) :: &
340 5 : "symdynmat", "symv1scf", "dvdb_add_lr", "interpolated"])
341 1 : NCF_CHECK(ncerr)
342 2 : NCF_CHECK(nctk_def_dpscalars(ncid, [character(len=nctk_slen) :: "qdamp"]))
343 :
344 : ! Define EPH arrays
345 : ncerr = nctk_def_arrays(ncid, [ &
346 : nctkarr_t('qpoint', "dp" , 'number_of_reduced_dimensions'), &
347 : nctkarr_t('emacro_cart', "dp", 'number_of_cartesian_directions, number_of_cartesian_directions'), &
348 : nctkarr_t('becs_cart', "dp", "number_of_cartesian_directions, number_of_cartesian_directions, number_of_atoms"), &
349 : nctkarr_t("eigenvalues_kq", "dp", "max_number_of_states, number_of_kpoints, number_of_spins"), &
350 : nctkarr_t('phfreqs', "dp", 'number_of_phonon_modes'), &
351 : nctkarr_t('phdispl_cart', "dp", 'complex, number_of_phonon_modes, number_of_phonon_modes'), &
352 : nctkarr_t('phdispl_red', "dp", 'complex, number_of_phonon_modes, number_of_phonon_modes'), &
353 : nctkarr_t("gkq_representation", "char", "character_string_length"), &
354 : nctkarr_t('gkq', "dp", &
355 : 'complex, max_number_of_states, max_number_of_states, number_of_phonon_modes, number_of_kpoints, number_of_spins') &
356 10 : ])
357 1 : NCF_CHECK(ncerr)
358 : ! Write data.
359 1 : NCF_CHECK(nctk_set_datamode(ncid))
360 : ncerr = nctk_write_iscalars(ncid, [character(len=nctk_slen) :: &
361 : "symdynmat", "symv1scf", "dvdb_add_lr", "interpolated"], &
362 9 : [dtset%symdynmat, dtset%symv1scf, dtset%dvdb_add_lr, interpolated])
363 1 : NCF_CHECK(ncerr)
364 2 : dvdb_qdamp = dvdb%qdamp
365 2 : NCF_CHECK(nctk_write_dpscalars(ncid, [character(len=nctk_slen) :: "qdamp"], dvdb_qdamp))
366 1 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "qpoint"), qpt))
367 1 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "emacro_cart"), dvdb%dielt))
368 1 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "becs_cart"), dvdb%zeff))
369 5 : ABI_MALLOC(eigens_kq, (ebands_kq%mband, nkpt, nsppol))
370 9 : do ik=1,nkpt
371 32 : kk = ebands_k%kptns(:,ik)
372 32 : kq = kk + qpt
373 : ! Find the index of the k+q point
374 8 : call findqg0(ikq, g0_k, kq, nkpt_kq, ebands_kq%kptns, [1,1,1])
375 57 : eigens_kq(:, ik, :) = ebands_kq%eig(:, ikq, :)
376 : end do
377 1 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "eigenvalues_kq"), eigens_kq))
378 1 : ABI_FREE(eigens_kq)
379 1 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "phfreqs"), phfrq))
380 1 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, 'phdispl_cart'), displ_cart))
381 1 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, 'phdispl_red'), displ_red))
382 1 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "gkq_representation"), "atom"))
383 : end if ! master
384 :
385 : else
386 0 : ABI_ERROR(sjoin("Invalid value for eph_task:", itoa(dtset%eph_task)))
387 : end if
388 :
389 : ! Loop over all 3*natom perturbations.
390 8 : do ipc=1,natom3
391 6 : idir = mod(ipc-1, 3) + 1
392 6 : ipert = (ipc - idir) / 3 + 1
393 6 : write(msg, '(a,2(i0,1x))') " Treating ipert, idir = ", ipert, idir
394 6 : call wrtout(std_out, msg, do_flush=.True.)
395 1371 : if (dtset%eph_task == 2) gkk = zero
396 :
397 12 : do spin=1,nsppol
398 1950 : if (dtset%eph_task == -2) gkq_atm = zero
399 :
400 : ! Set up local potential vlocal1 with proper dimensioning, from vtrial1 taking into account the spin.
401 : call rf_transgrid_and_pack(spin,nspden,psps%usepaw,cplex,nfftf,nfft,ngfft,gs_hamkq%nvloc,&
402 6 : pawfgr,mpi_enreg,dummy_vtrial,v1scf(:,:,:,ipc),vlocal,vlocal1(:,:,:,:,ipc))
403 :
404 : ! Continue to initialize the Hamiltonian
405 6 : call gs_hamkq%load_spin(spin,vlocal=vlocal,with_nonlocal=.true.)
406 :
407 : ! Allocate workspace for wavefunctions. Make npw larger than expected.
408 24 : ABI_MALLOC(bras, (2, mpw*nspinor, mband))
409 18 : ABI_MALLOC(kets, (2, mpw*nspinor, mband))
410 18 : ABI_MALLOC(h1_kets, (2, mpw*nspinor, mband))
411 :
412 : ! GKA: This little block used to be right after the perturbation loop
413 : ! Prepare application of the NL part.
414 6 : call rf_hamkq%init(cplex,gs_hamkq,ipert,has_e1kbsc=.true.)
415 6 : call rf_hamkq%load_spin(spin,vlocal1=vlocal1(:,:,:,:,ipc),with_nonlocal=.true.)
416 :
417 54 : do ik=1,nkpt
418 : ! Only do a subset a k-points
419 48 : if (.not. (ik >= my_kstart .and. ik <= my_kstop)) cycle
420 :
421 192 : kk = ebands_k%kptns(:,ik)
422 192 : kq = kk + qpt
423 : ! Find the index of the k+q point
424 48 : call findqg0(ikq, g0_k, kq, nkpt_kq, ebands_kq%kptns, [1,1,1])
425 :
426 : ! Copy u_k(G)
427 48 : istwf_k = wfd_k%istwfk(ik); npw_k = wfd_k%npwarr(ik)
428 48 : ABI_CHECK(mpw >= npw_k, "mpw < npw_k")
429 14664 : kg_k(:,1:npw_k) = wfd_k%kdata(ik)%kg_k
430 288 : do ib2=1,mband
431 288 : call wfd_k%copy_cg(ib2, ik, spin, kets(1,1,ib2))
432 : end do
433 :
434 : ! Copy u_kq(G)
435 48 : istwf_kq = wfd_kq%istwfk(ikq); npw_kq = wfd_kq%npwarr(ikq)
436 48 : ABI_CHECK(mpw >= npw_kq, "mpw < npw_kq")
437 13968 : kg_kq(:,1:npw_kq) = wfd_kq%kdata(ikq)%kg_k
438 288 : do ib1=1,mband_kq
439 288 : call wfd_kq%copy_cg(ib1, ikq, spin, bras(1,1,ib1))
440 : end do
441 :
442 : ! if PAW, one has to solve a generalized eigenproblem
443 : ! Be careful here because I will need sij_opt==-1
444 48 : gen_eigenpb = (psps%usepaw==1)
445 48 : sij_opt = 0; if (gen_eigenpb) sij_opt = 1
446 144 : ABI_MALLOC(gs1c, (2,npw_kq*nspinor*((sij_opt+1)/2)))
447 :
448 : ! GKA: Previous loop on 3*natom perturbations used to start here
449 : ! This call is not optimal because there are quantities in out that do not depend on idir,ipert
450 : call getgh1c_setup(gs_hamkq,rf_hamkq,dtset,psps,kk,kq,idir,ipert,& ! In
451 : cryst%natom,cryst%rmet,cryst%gprimd,cryst%gmet,istwf_k,& ! In
452 : npw_k,npw_kq,useylmgr1,kg_k,ylm_k,kg_kq,ylm_kq,ylmgr_kq_dum,& ! In
453 48 : dkinpw,nkpg,nkpg1,kpg_k,kpg_kq,kinpw_kq,ffnl_k,ffnl_kq,ph3d_k,ph3d_kq) ! Out
454 :
455 : ! Calculate dvscf * psi_k, results stored in h1_kets on the k+q sphere.
456 : ! Compute H(1) applied to GS wavefunction Psi(0)
457 288 : do ib2=1,mband
458 240 : eig0nk = ebands_k%eig(ib2,ik,spin)
459 : ! Use scissor shift on 0-order eigenvalue
460 240 : eshift = eig0nk - dtset%dfpt_sciss
461 :
462 : call getgh1c(berryopt0,kets(:,:,ib2),cwaveprj0,h1_kets(:,:,ib2),&
463 : grad_berry,gs1c,gs_hamkq,gvnlx1,idir,ipert, [eshift], mpi_enreg,1,optlocal,&
464 528 : optnl,opt_gvnlx1,rf_hamkq,sij_opt,tim_getgh1c,usevnl)
465 : end do
466 :
467 48 : ABI_FREE(kinpw_kq)
468 48 : ABI_FREE(kpg_k)
469 48 : ABI_FREE(kpg_kq)
470 48 : ABI_FREE(dkinpw)
471 48 : ABI_FREE(ffnl_k)
472 48 : ABI_FREE(ffnl_kq)
473 48 : ABI_FREE(gs1c)
474 48 : ABI_FREE(ph3d_k)
475 48 : ABI_SFREE(ph3d_kq)
476 :
477 : ! Calculate elphmat(j,i) = <psi_{k+q,j}|dvscf_q*psi_{k,i}> for this perturbation.
478 : ! The array eig1_k contains:
479 : !
480 : ! <u_(band,k+q)^(0)|H_(k+q,k)^(1)|u_(band,k)^(0)> (NC psps)
481 : ! <u_(band,k+q)^(0)|H_(k+q,k)^(1)-(eig0_k+eig0_k+q)/2.S^(1)|u_(band,k)^(0)> (PAW)
482 294 : do ib2=1,mband
483 1488 : do ib1=1,mband_kq
484 : call dotprod_g(dotr,doti,istwf_kq,npw_kq*nspinor,2,bras(1,1,ib1),h1_kets(1,1,ib2),&
485 1200 : mpi_enreg%me_g0,mpi_enreg%comm_spinorfft)
486 1200 : band = 2*ib2-1 + (spin-1) * 2 * mband
487 1440 : if (dtset%eph_task == 2) then
488 600 : gkk(band,ik,1,1,ib1) = dotr
489 600 : gkk(band+1,ik,1,1,ib1) = doti
490 : else
491 1800 : gkq_atm(:, ib1, ib2, ik) = [dotr, doti]
492 : end if
493 : end do ! ib1
494 : end do ! ib2
495 :
496 : end do ! ikpt
497 :
498 6 : ABI_FREE(bras)
499 6 : ABI_FREE(kets)
500 6 : ABI_FREE(h1_kets)
501 6 : call rf_hamkq%free()
502 :
503 12 : if (dtset%eph_task == -2) then
504 : ! Gather the k-points computed by all processes
505 3 : call xmpi_sum_master(gkq_atm, master, comm, ierr)
506 3 : if (i_am_master) then
507 : ! Write the netCDF file.
508 : ncerr = nf90_put_var(ncid, nctk_idname(ncid, "gkq"), gkq_atm, &
509 39 : start=[1, 1, 1, ipc, 1, spin], count=[2, mband, mband, 1, nkpt, 1])
510 3 : NCF_CHECK(ncerr)
511 : end if
512 : end if
513 :
514 : end do ! spin
515 :
516 8 : if (dtset%eph_task == 2) then
517 : ! Gather the k-points computed by all processes
518 3 : call xmpi_sum_master(gkk,master,comm,ierr)
519 : ! Init a gkk_t object
520 3 : call gkk_init(gkk,gkk2d,mband,nsppol,nkpt,1,1)
521 : ! Write the netCDF file.
522 3 : call appdig(ipc,dtfil%fnameabo_gkk,gkkfilnam)
523 3 : fname = strcat(gkkfilnam, ".nc")
524 3 : if (i_am_master) then
525 3 : NCF_CHECK_MSG(nctk_open_create(ncid, fname, xmpi_comm_self), "Creating GKK file")
526 3 : NCF_CHECK(cryst%ncwrite(ncid))
527 3 : NCF_CHECK(ebands_k%ncwrite(ncid))
528 3 : call gkk_ncwrite(gkk2d, qpt, 1.0_dp, ncid)
529 3 : NCF_CHECK(nf90_close(ncid))
530 : end if
531 : ! Free memory
532 3 : call gkk_free(gkk2d)
533 : end if
534 : end do ! ipc (loop over 3*natom atomic perturbations)
535 :
536 2 : call cwtime(cpu, wall, gflops, "stop")
537 2 : write(msg, '(2a)') " Computation of gkq matrix elements with ", trim(what)
538 2 : call wrtout(units, msg, do_flush=.True.)
539 2 : call wrtout(std_out, sjoin("cpu-time:", sec2str(cpu), ",wall-time:", sec2str(wall)), do_flush=.True.)
540 :
541 2 : if (dtset%eph_task == -2 .and. i_am_master) then
542 1 : NCF_CHECK(nf90_close(ncid))
543 : end if
544 :
545 : ! ===========
546 : ! Free memory
547 : ! ===========
548 2 : ABI_SFREE(gkk)
549 2 : ABI_SFREE(gkq_atm)
550 2 : ABI_FREE(displ_cart)
551 2 : ABI_FREE(displ_red)
552 2 : ABI_FREE(v1scf)
553 2 : ABI_FREE(vlocal1)
554 2 : ABI_FREE(gvnlx1)
555 2 : ABI_FREE(grad_berry)
556 2 : ABI_FREE(dummy_vtrial)
557 2 : ABI_FREE(ph1d)
558 2 : ABI_FREE(vlocal)
559 2 : ABI_FREE(kg_k)
560 2 : ABI_FREE(kg_kq)
561 2 : ABI_FREE(ylm_k)
562 2 : ABI_FREE(ylm_kq)
563 2 : call pawcprj_free(cwaveprj0)
564 2 : ABI_FREE(cwaveprj0)
565 :
566 2 : call gs_hamkq%free(); call wfd_k%free(); call wfd_kq%free()
567 :
568 10 : end subroutine eph_gkk
569 : !!***
570 :
571 : !----------------------------------------------------------------------
572 :
573 : !!****f* m_gkk/ncwrite_v1qnu
574 : !! NAME
575 : !! ncwrite_v1qnu
576 : !!
577 : !! FUNCTION
578 : !! Compute \delta V_{q,nu)(r) and dump results to netcdf file.
579 : !! This routine should be called by a single processor.
580 : !!
581 : !! INPUT
582 : !! dvdb<dbdb_type>=Database with the DFPT SCF potentials.
583 : !! dtset<dataset_type>= Input variables.
584 : !! ifc<ifc_type>=interatomic force constants and corresponding real space grid info.
585 : !! out_ncpath=Name of the netcdf file.
586 : !!
587 : !! OUTPUT
588 : !! Only writing
589 : !!
590 : !! SOURCE
591 :
592 0 : subroutine ncwrite_v1qnu(dvdb, dtset, ifc, out_ncpath)
593 :
594 : use m_bz_mesh, only : kpath_t
595 :
596 : !Arguments ------------------------------------
597 : class(dvdb_t),intent(inout) :: dvdb
598 : type(dataset_type),target,intent(in) :: dtset
599 : type(ifc_type),intent(in) :: ifc
600 : character(len=*),intent(in) :: out_ncpath
601 :
602 : !Local variables-------------------------------
603 : !scalars
604 : integer,parameter :: master = 0, qptopt1 = 1
605 : integer :: db_iqpt, cplex, nfft, comm, ip, idir, ipert, my_rank, interpolated, comm_rpt, ncid, ncerr
606 : integer :: iq, nu, iatom, ii, jj, kk
607 : real(dp) :: inv_qepsq, qtau, phre, phim, rtmp
608 : logical :: with_lr_model
609 0 : type(kpath_t) :: qpath
610 : !arrays
611 : integer :: ngfft(18), units(2)
612 0 : real(dp) :: phfreqs(dvdb%natom3),qpt(3)
613 0 : real(dp) :: displ_cart(2,3, dvdb%cryst%natom, dvdb%natom3), displ_red(2,dvdb%natom3,dvdb%natom3)
614 0 : real(dp),allocatable :: v1scf(:,:,:,:), v1_qnu(:,:,:,:), v1lr_atm(:,:,:,:), v1lr_qnu(:,:,:,:)
615 0 : real(dp) :: bounds(3,6), qpt_red(3), qpt_cart(3), glr(3), values(dvdb%natom3)
616 :
617 : !************************************************************************
618 :
619 : ! +0.50000 +0.50000 +0.50000 # L
620 : ! +0.00000 +0.00000 +0.00000 # $\Gamma$
621 : ! +0.50000 +0.00000 +0.50000 # X
622 : ! +0.50000 +0.25000 +0.75000 # W
623 : ! +0.37500 +0.37500 +0.75000 # K
624 : ! +0.00000 +0.00000 +0.00000 # $\Gamma$
625 : ! +0.37500 +0.37500 +0.75000 # K
626 :
627 : ! +0.62500 +0.25000 +0.62500 # U
628 : ! +0.50000 +0.50000 +0.50000 # L
629 : ! +0.37500 +0.37500 +0.75000 # K
630 : ! +0.62500 +0.25000 +0.62500 # U
631 : ! +0.50000 +0.00000 +0.50000 # X
632 :
633 0 : bounds(:, 1) = tol3 * [+0.50000, +0.50000, +0.50000] ! # L
634 0 : bounds(:, 2) = tol3 * [+0.00000, +0.00000, +0.00000] ! # $\Gamma$
635 0 : bounds(:, 3) = tol3 * [+0.50000, +0.00000, +0.50000] ! # X
636 0 : bounds(:, 4) = tol3 * [+0.37500, +0.37500, +0.75000] ! # K
637 0 : bounds(:, 5) = tol3 * [+0.00000, +0.00000, +0.00000] ! # $\Gamma$
638 0 : bounds(:, 6) = tol3 * [+0.50000, +0.25000, +0.75000] ! # W
639 :
640 0 : call qpath%init(bounds, dvdb%cryst%gprimd, dtset%ndivsm)
641 :
642 0 : units = [std_out, ab_out]
643 :
644 0 : do iq=1,qpath%npts
645 0 : qpt_red = qpath%points(:, iq)
646 0 : qpt_cart = two_pi * matmul(dvdb%cryst%gprimd, qpt_red)
647 0 : inv_qepsq = one / dot_product(qpt_cart, matmul(ifc%dielt, qpt_cart))
648 0 : call ifc%fourq(dvdb%cryst, qpt_red, phfreqs, displ_cart)
649 0 : do nu=1, dvdb%natom3
650 0 : glr = zero
651 0 : do iatom=1, dvdb%cryst%natom
652 : ! Phase factor exp(-i (q+G) . tau)
653 0 : qtau = - two_pi * dot_product(qpt_red, dvdb%cryst%xred(:,iatom))
654 0 : phre = cos(qtau); phim = sin(qtau)
655 0 : do jj=1,3
656 0 : do ii=1,3
657 0 : do kk=1,3
658 0 : rtmp = dvdb%qstar(ii, jj, kk, iatom) * qpt_cart(ii) * qpt_cart(jj)
659 0 : glr(1) = glr(1) + rtmp * (displ_cart(1, kk, iatom, nu) * phre - displ_cart(2, kk, iatom, nu) * phim)
660 0 : glr(2) = glr(2) + rtmp * (displ_cart(2, kk, iatom, nu) * phre + displ_cart(1, kk, iatom, nu) * phre)
661 : end do
662 : end do
663 : end do
664 : end do
665 0 : glr = half * (glr / inv_qepsq) * (four_pi / dvdb%cryst%ucvol)
666 0 : values(nu) = (glr(1) ** 2 + glr(2) ** 2) / (two * phfreqs(nu))
667 : end do ! nu
668 0 : write(std_out, "(i0, 4(f9.6), /, (es18.6, 1x))") iq, qpt_red, phfreqs(nu), (values(nu), nu=1, 3*dvdb%natom)
669 : end do ! iqpt
670 :
671 0 : call qpath%free()
672 : return
673 :
674 : my_rank = xmpi_comm_rank(dvdb%comm)
675 : comm = dvdb%comm
676 : qpt = dtset%qptn
677 :
678 : call wrtout(std_out, sjoin(" Writing Delta V_{q,nu)(r) potentials to file:", out_ncpath), do_flush=.True.)
679 : call wrtout(units, sjoin(ch10, "- Results stored in: ", out_ncpath))
680 : call wrtout(std_out, sjoin(" Using qpt:", ktoa(qpt)))
681 : !call wrtout(units, " Use `abiopen.py out_V1QAVG.nc -e` to visualize results")
682 : call dvdb%print([std_out], "", 0)
683 :
684 : ! Define FFT mesh
685 : ngfft = dvdb%ngfft
686 : nfft = product(ngfft(1:3))
687 :
688 : if (dtset%eph_task == -16) then
689 : call wrtout(units, " Assuming q-point already in the DVDB file. No interpolation.")
690 : interpolated = 0
691 :
692 : else if (dtset%eph_task == +16) then
693 : call wrtout(units, " Using Fourier interpolation.")
694 : comm_rpt = xmpi_comm_self
695 : call dvdb%ftinterp_setup(dtset%ddb_ngqpt, qptopt1, 1, dtset%ddb_shiftq, nfft, ngfft, comm_rpt)
696 : interpolated = 1
697 : else
698 : ABI_ERROR(sjoin("Invalid value for eph_task:", itoa(dtset%eph_task)))
699 : end if
700 :
701 : with_lr_model = .True.
702 :
703 : ! Create netcdf file.
704 : if (my_rank == master) then
705 : NCF_CHECK(nctk_open_create(ncid, out_ncpath, comm))
706 : NCF_CHECK(dvdb%cryst%ncwrite(ncid))
707 :
708 : ! Add other dimensions.
709 : ncerr = nctk_def_dims(ncid, [ &
710 : nctkdim_t("nfft", nfft), nctkdim_t("nspden", dvdb%nspden), &
711 : nctkdim_t("natom3", 3 * dvdb%cryst%natom)], defmode=.True.)
712 : NCF_CHECK(ncerr)
713 :
714 : ! Define arrays
715 : ncerr = nctk_def_arrays(ncid, [ &
716 : nctkarr_t("ngfft", "int", "three"), &
717 : nctkarr_t("qpt", "dp", "three"), &
718 : nctkarr_t("phfreqs", "dp", "natom3"), &
719 : nctkarr_t("displ_cart", "dp", "two, natom3, natom3"), &
720 : nctkarr_t("v1_qnu", "dp", "two, nfft, nspden, natom3")])
721 : NCF_CHECK(ncerr)
722 :
723 : if (with_lr_model) then
724 : NCF_CHECK(nctk_def_arrays(ncid, [nctkarr_t("v1lr_qnu", "dp", "two, nfft, nspden, natom3")]))
725 : end if
726 :
727 : NCF_CHECK(nctk_set_datamode(ncid))
728 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "ngfft"), ngfft(1:3)))
729 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "qpt"), qpt))
730 : end if
731 :
732 : ABI_MALLOC(v1_qnu, (2, nfft, dvdb%nspden, dvdb%natom3))
733 : if (with_lr_model) then
734 : ABI_MALLOC(v1lr_atm, (2, nfft, dvdb%nspden, dvdb%natom3))
735 : ABI_MALLOC(v1lr_qnu, (2, nfft, dvdb%nspden, dvdb%natom3))
736 : end if
737 :
738 : ! Get phonon freqs and displacemented for this q-point.
739 : call ifc%fourq(dvdb%cryst, qpt, phfreqs, displ_cart, out_displ_red=displ_red)
740 :
741 : if (interpolated == 0) then
742 : ! Find the index of the q-point in the DVDB.
743 : db_iqpt = dvdb%findq(qpt)
744 : if (db_iqpt /= -1) then
745 : ! Read or reconstruct the dvscf potentials for all 3*natom perturbations.
746 : ! This call allocates v1scf(cplex, nfft, nspden, 3*natom))
747 : call dvdb%readsym_allv1(db_iqpt, cplex, nfft, ngfft, v1scf, comm)
748 : else
749 : ABI_ERROR(sjoin("Cannot find q-point:", ktoa(qpt), "in DVDB file"))
750 : end if
751 : else
752 :
753 : cplex = 2
754 : ABI_MALLOC(v1scf, (cplex, nfft, dvdb%nspden, dvdb%my_npert))
755 : call dvdb%ftinterp_qpt(qpt, nfft, ngfft, v1scf, dvdb%comm_rpt)
756 : end if
757 :
758 : ! Compute scattering potential the in phonon representations instead of atomic one.
759 : ! v1_qnu = \sum_{ka} phdispl{ka}(q,nu) D_{ka,q} V_scf(r)
760 : ! NOTE: prefactor 1/sqrt(2 w(q,nu)) is not included in the potentials saved to file.
761 : ! v1_qnu(2, nfft, nspden, natom3), v1scf(cplex, nfft, nspden, natom3)
762 : call ephtk_v1atm_to_vqnu(cplex, nfft, dvdb%nspden, dvdb%natom3, v1scf, displ_red, v1_qnu)
763 :
764 : if (with_lr_model) then
765 : ! Compute LR model in the atomic representation then compute phonon representation in v1lr_qnu.
766 : v1lr_atm = zero
767 : do idir=1,3
768 : do ipert=1,dvdb%natom
769 : ip = (ipert - 1) * 3 + idir
770 : call dvdb%get_v1r_long_range(qpt, idir, ipert, nfft, ngfft, v1lr_atm(:,:,1,ip))
771 : if (dvdb%nspden == 2) v1lr_atm(:,:,2,ip) = v1lr_atm(:,:,1,ip)
772 : end do
773 : end do
774 : call ephtk_v1atm_to_vqnu(2, nfft, dvdb%nspden, dvdb%natom3, v1lr_atm, displ_red, v1lr_qnu)
775 : end if
776 :
777 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "phfreqs"), phfreqs))
778 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "displ_cart"), displ_cart))
779 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "v1_qnu"), v1_qnu))
780 : if (with_lr_model) then
781 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "v1lr_qnu"), v1lr_qnu))
782 : end if
783 :
784 : ABI_FREE(v1scf)
785 : ABI_FREE(v1_qnu)
786 : ABI_SFREE(v1lr_atm)
787 : ABI_SFREE(v1lr_qnu)
788 :
789 : NCF_CHECK(nf90_close(ncid))
790 : call dvdb%close()
791 :
792 : call wrtout(std_out, "dvqnu file written", do_flush=.True.)
793 :
794 0 : end subroutine ncwrite_v1qnu
795 : !!***
796 :
797 : !----------------------------------------------------------------------
798 :
799 : !!****f* m_gkk/find_mpw
800 : !! NAME
801 : !! find_mpw
802 : !!
803 : !! FUNCTION
804 : !! Look at all k-points and spins to find the maximum number of plane waves.
805 : !!
806 : !! INPUTS
807 : !!
808 : !! OUTPUT
809 : !!
810 : !! SOURCE
811 :
812 4 : subroutine find_mpw(mpw, kpts, nsppol, nkpt, gmet, ecut, comm)
813 :
814 : !Arguments ------------------------------------
815 : !scalars
816 : integer,intent(out) :: mpw
817 : integer,intent(in) :: nsppol, nkpt, comm
818 : real(dp),intent(in) :: ecut
819 : !arrays
820 : real(dp),intent(in) :: kpts(3,nkpt), gmet(3,3)
821 :
822 : !Local variables ------------------------------
823 : !scalars
824 : integer :: my_rank, cnt, nproc, ierr, ispin, ikpt, my_mpw, onpw
825 4 : integer,allocatable :: gtmp(:,:)
826 : real(dp) :: kpt(3)
827 :
828 : !************************************************************************
829 :
830 4 : my_rank = xmpi_comm_rank(comm); nproc = xmpi_comm_size(comm)
831 :
832 4 : mpw = 0; cnt=0
833 8 : do ispin=1,nsppol
834 40 : do ikpt=1,nkpt
835 32 : cnt = cnt + 1; if (mod(cnt, nproc) /= my_rank) cycle
836 128 : kpt = kpts(:,ikpt)
837 32 : call get_kg(kpt,1,ecut,gmet,onpw,gtmp)
838 32 : ABI_FREE(gtmp)
839 36 : mpw = max(mpw, onpw)
840 : end do
841 : end do
842 4 : my_mpw = mpw; call xmpi_max(my_mpw, mpw, comm, ierr)
843 :
844 4 : end subroutine find_mpw
845 : !!***
846 :
847 : end module m_gkk
848 : !!***
|
