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