Branch data Line data Source code
1 : : !!****m* ABINIT/m_ddk
2 : : !! NAME
3 : : !! m_ddk
4 : : !!
5 : : !! FUNCTION
6 : : !! Objects and methods to extract data from DDK files.
7 : : !! The DDK files are binary (soon also netcdf) files with Hamiltonian derivatives
8 : : !! wrt k, and the corresponding wave functions
9 : : !!
10 : : !! COPYRIGHT
11 : : !! Copyright (C) 2016-2026 ABINIT group (MG, HM, MJV)
12 : : !! This file is distributed under the terms of the
13 : : !! GNU General Public License, see ~abinit/COPYING
14 : : !! or http://www.gnu.org/copyleft/gpl.txt .
15 : : !! For the initials of contributors, see ~abinit/doc/developers/contributors.txt.
16 : : !!
17 : : !! SOURCE
18 : :
19 : : #if defined HAVE_CONFIG_H
20 : : #include "config.h"
21 : : #endif
22 : :
23 : : #include "abi_common.h"
24 : :
25 : : module m_ddk
26 : :
27 : : use defs_basis
28 : : use m_abicore
29 : : use m_errors
30 : : use m_xmpi
31 : : use m_nctk
32 : : use m_hdr
33 : : use m_dtset
34 : : use m_krank
35 : : use m_crystal
36 : : use m_mpinfo
37 : : use m_cgtools
38 : : use m_hamiltonian
39 : : use m_initylmg
40 : : use m_pawcprj
41 : : use m_getgh1c
42 : : use netcdf
43 : :
44 : : use m_fstrings, only : strcat, sjoin, itoa, ktoa
45 : : use m_io_tools, only : iomode_from_fname
46 : : use m_time, only : cwtime, cwtime_report
47 : : use defs_abitypes, only : MPI_type
48 : : use defs_datatypes, only : pseudopotential_type
49 : : use m_vkbr, only : vkbr_t, nc_ihr_comm, vkbr_init, vkbr_free
50 : : use m_pawtab, only : pawtab_type
51 : : use m_wfk, only : wfk_read_ebands !, wfk_read_h1mat
52 : : use m_wfd, only : wfd_t, wave_t
53 : : use m_ebands, only : ebands_t
54 : :
55 : : implicit none
56 : :
57 : : private
58 : : !!***
59 : :
60 : : public :: ddk_red2car ! Convert band velocities from cartesian to reduced coordinates
61 : :
62 : : !!***
63 : :
64 : : type, private :: ham_targets_t
65 : : real(dp),allocatable :: ffnlk(:,:,:,:), ffnl1(:,:,:,:)
66 : : real(dp),allocatable :: kpg_k(:,:), kpg1_k(:,:)
67 : : real(dp),allocatable :: ph3d(:,:,:), ph3d1(:,:,:)
68 : : real(dp),allocatable :: dkinpw(:), kinpw1(:)
69 : : contains
70 : : procedure :: free => ham_targets_free ! Free memory.
71 : : end type ham_targets_t
72 : :
73 : :
74 : : !!****t* m_ddk/ddkop_t
75 : : !! NAME
76 : : !! ddkop_t
77 : : !!
78 : : !! FUNCTION
79 : : !! This object provides a simplified interface to compute matrix elements of the
80 : : !! velocity operator with the DFPT routines.
81 : : !!
82 : : !! SOURCE
83 : :
84 : : type,public :: ddkop_t
85 : :
86 : : integer :: ipert
87 : : ! Perturbation type: natom + 1
88 : :
89 : : integer :: inclvkb
90 : : ! Option for calculating the matrix elements of [Vnl,r].
91 : : ! 0 to exclude commutator, 2 to include it
92 : :
93 : : integer :: usepaw
94 : : ! 0 for NC, 1 for PAW
95 : :
96 : : integer :: mpw
97 : : ! Maximum number of plane-waves over k-points (used to dimension arrays)
98 : :
99 : : real(dp) :: kpoint(3)
100 : : ! K-point (set in setup_spin_kpoint)
101 : :
102 : : real(dp) :: eig0nk
103 : :
104 : : real(dp) :: dfpt_sciss = zero
105 : :
106 : : real(dp) :: rprimd(3,3)
107 : :
108 : : type(MPI_type),pointer :: mpi_enreg => null()
109 : :
110 : : type(gs_hamiltonian_type) :: gs_hamkq(3)
111 : :
112 : : type(rf_hamiltonian_type) :: rf_hamkq(3)
113 : :
114 : : type(ham_targets_t), private :: htg(3)
115 : : ! Store arrays targetted by the hamiltonians.
116 : :
117 : : real(dp), allocatable :: gh1c(:,:,:)
118 : : !gh1c, (2, npw_k*nspinor, 3))
119 : :
120 : : real(dp), allocatable :: gs1c(:,:,:)
121 : : ! gs1c, (2, npw_k*nspinor, 3*psps%usepaw))
122 : :
123 : : contains
124 : :
125 : : procedure :: init => ddkop_init
126 : : ! Build object
127 : :
128 : : procedure :: setup_spin_kpoint => ddkop_setup_spin_kpoint
129 : : ! Prepare application of dH/dk for given spin, k-point.
130 : :
131 : : procedure :: apply => ddkop_apply
132 : : ! Apply dH/dk to input wavefunction.
133 : :
134 : : procedure :: get_braket => ddkop_get_braket
135 : : ! Compute matrix element (complex results) in cartesian coords.
136 : :
137 : : procedure :: get_ihr_comm => ddkop_get_ihr_comm
138 : : ! Compute matrix elements of i[H,r] (complex results) in cartesian coords.
139 : :
140 : : procedure :: get_vdiag => ddkop_get_vdiag
141 : : ! Compute diagonal matrix element (real) in cartesian coords.
142 : :
143 : : procedure :: get_vnondiag => ddkop_get_vnondiag
144 : : ! Compute off diagonal matrix elements in cartesian coords.
145 : :
146 : : procedure :: free => ddkop_free
147 : : ! Free memory.
148 : :
149 : : end type ddkop_t
150 : : !!***
151 : :
152 : : !!****t* m_ddk/ddkstore_t
153 : : !! NAME
154 : : !! ddkstore_t
155 : : !!
156 : : !! FUNCTION
157 : : !! This object stores the matrix elements of the velocity operator computed with the DFPT routines.
158 : : !!
159 : : !! SOURCE
160 : :
161 : : type,public :: ddkstore_t
162 : :
163 : : integer :: bmin = 1, bmax = -1
164 : : ! Min and max band index
165 : :
166 : : character(len=50) :: mode = "reduced"
167 : : ! "cart" or "reduced"
168 : :
169 : : logical :: only_diago = .False.
170 : : ! True if we are computing only the diagonal elements
171 : :
172 : : real(dp),allocatable :: dipoles(:,:,:,:,:,:)
173 : : ! (3, 2, mband, mband, nkpt, nsppol))
174 : :
175 : : real(dp),allocatable :: vdiago(:,:,:,:)
176 : : ! (3, bmin:bmax, nkpt, nsppol)
177 : :
178 : : real(dp),allocatable :: vmat(:,:,:,:,:,:)
179 : : ! (2, 3, bmin:bmax, bmin:bmax, nkpt, nsppol))
180 : :
181 : : contains
182 : :
183 : : procedure :: compute_ddk => ddkstore_compute_ddk
184 : : ! Calculate DDK matrix elements (diago or full b,b' matrix).
185 : : ! Return results in datatype. Optionally, save results to disk in EVK format.
186 : :
187 : : procedure :: free => ddkstore_free
188 : : ! Free memory.
189 : :
190 : : end type ddkstore_t
191 : : !!***
192 : :
193 : : CONTAINS
194 : :
195 : : !----------------------------------------------------------------------
196 : :
197 : : !!****f* m_ddk/ddkstore_compute_ddk
198 : : !! NAME
199 : : !! ddkstore_compute_ddk
200 : : !!
201 : : !! FUNCTION
202 : : !! Calculate the DDK matrix elements using the commutator formulation.
203 : : !!
204 : : !! INPUTS
205 : : !! prefix: Prefix for output EVK file. Empty if output files are not wanted
206 : : !!
207 : : !! SOURCE
208 : :
209 : 5 : subroutine ddkstore_compute_ddk(ds, wfk_path, prefix, dtset, psps, pawtab, ngfftc, comm)
210 : :
211 : : !Arguments ------------------------------------
212 : : !scalars
213 : : class(ddkstore_t),intent(inout) :: ds
214 : : character(len=*),intent(in) :: wfk_path, prefix
215 : : integer,intent(in) :: comm
216 : : type(dataset_type),intent(in) :: dtset
217 : : type(pseudopotential_type),intent(in) :: psps
218 : : type(pawtab_type),intent(in) :: pawtab(psps%ntypat*psps%usepaw)
219 : : !arrays
220 : : integer,intent(in) :: ngfftc(18)
221 : :
222 : : !Local variables ------------------------------
223 : : !scalars
224 : : integer,parameter :: master = 0
225 : : integer :: mband, nbcalc, nsppol, ib_v, ib_c, mpw, spin, nspinor, nkpt, nband_k, npw_k
226 : : integer :: ii, ik, bmin, bmax, istwf_k, idir, my_rank, nproc, ierr, bstop
227 : : real(dp) :: cpu, wall, gflops, cpu_all, wall_all, gflops_all
228 : : integer :: ncerr, ncid
229 : : character(len=500) :: msg
230 : : character(len=fnlen) :: fname
231 : : logical :: write_ncfile
232 : 5 : type(wfd_t) :: wfd
233 : 5 : type(vkbr_t) :: vkbr
234 : 5 : type(ebands_t) :: ebands
235 [ + + + + : 260 : type(crystal_t) :: cryst
+ + + + +
+ + + + +
+ + + + ]
236 : 5 : type(hdr_type) :: tmp_hdr, hdr
237 [ + + - + : 50 : type(ddkop_t) :: ddkop
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ + + - +
- + + + -
+ - + - +
- + - + -
+ - + -
+ ]
238 : : type(wave_t),pointer :: wave_v, wave_c
239 : : !arrays
240 : 5 : integer,allocatable :: distrib_mat(:,:,:,:), distrib_diago(:,:,:), nband(:,:), kg_k(:,:)
241 : 10 : logical,allocatable :: bks_mask(:,:,:), keep_ur(:,:,:)
242 : : real(dp) :: kpt(3), vv(2, 3)
243 : 5 : real(dp),allocatable :: cg_c(:,:), cg_v(:,:)
244 : : complex(dp) :: vg(3), vr(3)
245 : 5 : complex(gwp),allocatable :: ihrc(:,:), ug_c(:), ug_v(:)
246 : 5 : type(pawcprj_type),allocatable :: cwaveprj(:,:)
247 : : !************************************************************************
248 : :
249 : 5 : my_rank = xmpi_comm_rank(comm); nproc = xmpi_comm_size(comm)
250 : :
251 [ + - + + ]: 15 : if (my_rank == master) call wrtout([std_out, ab_out], " Computation of velocity matrix elements (DDK)", newlines=1)
252 : :
253 [ - + ]: 5 : ABI_CHECK(psps%usepaw == 0, "PAW not implemented")
254 : :
255 : : ! Get ebands and hdr from WFK file.
256 [ - + - + : 5 : ebands = wfk_read_ebands(wfk_path, comm, out_hdr=hdr)
- + - + -
+ - + - +
- + - + -
+ - + ]
257 [ - + - + : 5 : cryst = hdr%get_crystal()
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + ]
258 : :
259 : : ! Extract important dimensions from hdr%
260 : 5 : nkpt = hdr%nkpt
261 : 5 : nsppol = hdr%nsppol
262 : 5 : nspinor = hdr%nspinor
263 : 5 : mband = hdr%mband
264 : :
265 : : ! Define band range
266 : : ! TODO: Perhaps one should allocate output arrays using bmin:bmax
267 : : ! and allow for nc output only if bmin == 1 and bmax == mband
268 [ + - ]: 5 : if (ds%bmax == -1) ds%bmax = mband
269 : :
270 [ + - + - : 5 : if (ds%bmin < 1 .or. ds%bmin > mband .or. ds%bmax > mband .or. ds%bmin > ds%bmax) then
+ - - + ]
271 : 0 : ABI_ERROR(sjoin("Invalid value for bmin, bmax", itoa(ds%bmin), itoa(ds%bmax), "with mband:", itoa(mband)))
272 : : end if
273 : :
274 : 5 : bmin = ds%bmin; bmax = ds%bmax
275 : 5 : nbcalc = bmax - bmin + 1
276 : 5 : write_ncfile = len_trim(prefix) > 0
277 [ - + + - : 5 : if (write_ncfile .and. .not. (bmin == 1 .and. bmax == mband) ) then
- + ]
278 : 0 : write_ncfile = .False.
279 : 0 : ABI_WARNING("Cannot write ncfile if .not. (bmin == 1 .and. bmax == mband)")
280 : : end if
281 : :
282 [ + - ]: 5 : if (my_rank == master) then
283 : 5 : write(ab_out, "(a)")" Parameters extracted from the Abinit header:"
284 : 5 : write(ab_out, "(a, f5.1)") ' ecut: ', hdr%ecut
285 : 5 : write(ab_out, "(a, i0)") ' nkpt: ', nkpt
286 : 5 : write(ab_out, "(a, i0)") ' mband: ', mband
287 : 5 : write(ab_out, "(a, i0)") ' nsppol: ', nsppol
288 : 5 : write(ab_out, "(a, i0)") ' nspinor: ', nspinor
289 : 5 : write(ab_out, "(a, i0)") ' useylm: ', dtset%useylm
290 : 5 : write(ab_out, "(a, i0)") ' inclvkb: ', dtset%inclvkb
291 : 5 : write(ab_out, "(2(a, i0))")' bmin: ', bmin, ", bmax: ", bmax
292 [ + + ]: 5 : if (ds%only_diago) then
293 : 1 : write(ab_out, "(a)")' Computing diagonal matrix elements only'
294 : : else
295 : 4 : write(ab_out, "(a)")' Computing diagonal and off-diagonal matrix elements'
296 : : end if
297 : 5 : write(ab_out, "(2(a, i0))")' Between band index bmin: ', bmin, ", bmax: ", bmax
298 : 5 : write(ab_out, "(a)")""
299 : : end if
300 : :
301 : : ! Create distribution of the wavefunctions mask.
302 [ - + + - : 20 : ABI_MALLOC(nband, (nkpt, nsppol))
- + + - +
- + - - +
- + ]
303 [ - + + - : 25 : ABI_MALLOC(keep_ur, (mband, nkpt, nsppol))
- + + - -
+ + - + -
+ - - + -
+ ]
304 [ - + + - : 20 : ABI_MALLOC(bks_mask, (mband, nkpt, nsppol))
- + + - -
+ + - - +
- + ]
305 [ + + + + : 2078 : keep_ur = .false.; bks_mask = .false.; nband = mband
+ + + + +
+ + + + +
+ + ]
306 : :
307 [ + + ]: 5 : if (ds%only_diago) then
308 : : ! Distribute k-points, spin and (b, b) diagonal over MPI processors.
309 [ - + + - : 5 : ABI_MALLOC(distrib_diago, (bmin:bmax, nkpt, nsppol))
- + + - -
+ + - + -
+ - - + -
+ ]
310 [ + + + + : 74 : distrib_diago = -1
+ + ]
311 : :
312 : : ! Create bks_mask to load the wavefunctions.
313 : 1 : ii = 0
314 [ + + ]: 2 : do spin=1,nsppol
315 [ + + ]: 10 : do ik=1,nkpt
316 [ + + ]: 73 : do ib_v=bmin,bmax
317 [ - + ]: 64 : ii = ii + 1; if (mod(ii, nproc) /= my_rank) cycle ! MPI parallelism.
318 : 64 : distrib_diago(ib_v, ik, spin) = my_rank
319 : 72 : bks_mask(ib_v, ik, spin) = .true.
320 : : end do
321 : : end do
322 : : end do
323 [ + + + + : 74 : call wrtout(std_out, sjoin(" Rank: ", itoa(my_rank), "will treat", itoa(count(distrib_diago == my_rank))))
+ + + - ]
324 : :
325 : : else
326 : : ! Distribute k-points, spin and (b, b') pairs over the processors
327 [ - + + - : 24 : ABI_MALLOC(distrib_mat, (bmin:bmax, bmin:bmax, nkpt, nsppol))
- + + - -
+ + - - +
+ - + - +
- - + - +
- + ]
328 : 4 : call xmpi_distab(nproc, distrib_mat)
329 : :
330 : : ! Create bks_mask to load the wavefunctions
331 [ + + ]: 8 : do spin=1,nsppol
332 [ + + ]: 102 : do ik=1,nkpt
333 : : ! Loop over v bands
334 [ + + ]: 910 : do ib_v=bmin,bmax
335 : : ! Loop over c bands
336 [ + + ]: 8002 : do ib_c=bmin,bmax
337 [ + - ]: 7908 : if (distrib_mat(ib_c, ib_v, ik, spin) == my_rank) then
338 : 7096 : bks_mask(ib_v, ik, spin) = .true.
339 : 7096 : bks_mask(ib_c, ik, spin) = .true.
340 : : end if
341 : : end do
342 : : end do
343 : : end do
344 : : end do
345 : :
346 [ + + + + : 8010 : call wrtout(std_out, sjoin(" Rank: ", itoa(my_rank), "will treat", itoa(count(distrib_mat == my_rank))))
+ + + + +
- ]
347 : : end if
348 : :
349 : : ! Initialize distributed wavefunctions object
350 : : call wfd%init(cryst, pawtab, psps, keep_ur, mband, nband, nkpt, nsppol,&
351 : : bks_mask, dtset%nspden, nspinor, hdr%ecut, dtset%ecutsm, dtset%dilatmx, ebands%istwfk, ebands%kptns,&
352 : 5 : ngfftc, dtset%nloalg, dtset%prtvol, dtset%pawprtvol, comm)
353 : :
354 [ - + ]: 5 : ABI_FREE(bks_mask)
355 [ - + ]: 5 : ABI_FREE(keep_ur)
356 [ - + ]: 5 : ABI_FREE(nband)
357 : :
358 [ + + ]: 10 : call wfd%print([std_out], header="Wavefunctions on the k-points grid")
359 : :
360 : : ! Read wavefunctions from WFK file.
361 : 5 : call wfd%read_wfk(wfk_path, iomode_from_fname(wfk_path))
362 : :
363 : : ! Allocate workspace arrays
364 [ + + ]: 107 : mpw = maxval(wfd%npwarr)
365 [ - + + - : 15 : ABI_MALLOC(kg_k, (3, mpw))
+ - + - -
+ - + ]
366 [ - + + - : 15 : ABI_MALLOC(ug_c, (mpw*nspinor))
+ - + - -
+ - + ]
367 [ - + + - : 10 : ABI_MALLOC(ug_v, (mpw*nspinor))
- + - + ]
368 [ + - ]: 5 : if (dtset%useria /= 666) then
369 [ - + + - : 15 : ABI_MALLOC(cg_c, (2, mpw*nspinor))
+ - + - -
+ - + ]
370 [ - + + - : 10 : ABI_MALLOC(cg_v, (2, mpw*nspinor))
- + - + ]
371 : : end if
372 : :
373 [ - + - + ]: 5 : ABI_MALLOC(cwaveprj, (0, 0))
374 [ - + + - : 69485 : ABI_CALLOC(ds%dipoles, (3, 2, bmin:bmax, bmin:bmax, nkpt, nsppol))
- + + - -
+ + - - +
+ - + - +
- - + - +
- + + + +
+ + + + +
+ + + + ]
375 [ - + + - : 15 : ABI_MALLOC(ihrc, (3, nspinor**2))
+ - + - -
+ - + -
+ ]
376 : :
377 [ + + ]: 5 : if (ds%only_diago) then
378 [ - + + - : 270 : ABI_CALLOC(ds%vdiago, (3, bmin:bmax, nkpt, nsppol))
- + + - -
+ + - + -
+ - - + -
+ - + + +
+ + + + +
+ ]
379 : : else
380 [ - + + - : 71890 : ABI_CALLOC(ds%vmat, (2, 3, bmin:bmax, bmin:bmax, nkpt, nsppol))
- + + - -
+ + - - +
+ - - + -
+ - + + +
+ + + + +
+ + + +
+ ]
381 : : end if
382 : :
383 [ + - ]: 5 : if (dtset%useria /= 666) then
384 : 5 : call ddkop%init(dtset, cryst, pawtab, psps, wfd%mpi_enreg, mpw, wfd%ngfft)
385 : : !if (my_rank == master) call ddkop%print(ab_out)
386 : : end if
387 : :
388 : 5 : call cwtime(cpu_all, wall_all, gflops_all, "start")
389 : :
390 [ + + ]: 10 : do spin=1,nsppol
391 [ + + ]: 112 : do ik=1,nkpt
392 : :
393 : : ! Only do a subset a k-points
394 [ + + ]: 102 : if (ds%only_diago) then
395 [ + - - + : 8 : if (all(distrib_diago(:, ik, spin) /= my_rank)) cycle
+ - ]
396 : : else
397 [ + - + - : 94 : if (all(distrib_mat(bmin:bmax, bmin:bmax, ik, spin) /= my_rank)) cycle
- + + - ]
398 : : end if
399 : 102 : call cwtime(cpu, wall, gflops, "start")
400 : :
401 : 102 : nband_k = wfd%nband(ik, spin)
402 : 102 : istwf_k = wfd%istwfk(ik)
403 : 102 : npw_k = wfd%npwarr(ik)
404 [ + + ]: 408 : kpt = wfd%kibz(:,ik)
405 [ + + + + ]: 220558 : kg_k(:,1:npw_k) = wfd%kdata(ik)%kg_k
406 : :
407 [ + - ]: 102 : if (dtset%useria /= 666) then
408 : 102 : call ddkop%setup_spin_kpoint(dtset, cryst, psps, spin, kpt, istwf_k, npw_k, kg_k)
409 : : else
410 : : ! Allocate KB form factors
411 : : ! Prepare term i <n,k|[Vnl,r]|n"k>
412 [ # # ]: 0 : if (dtset%inclvkb /= 0) call vkbr_init(vkbr, cryst, psps, dtset%inclvkb, istwf_k, npw_k, kpt, kg_k)
413 : : end if
414 : :
415 : : ! Loop over bands
416 [ + + ]: 978 : do ib_v=bmin,bmax
417 [ + + ]: 876 : if (ds%only_diago) then
418 [ + - ]: 64 : if (distrib_diago(ib_v,ik,spin) /= my_rank) cycle
419 : : else
420 [ + - - + : 812 : if (all(distrib_mat(:,ib_v,ik,spin) /= my_rank)) cycle
- + ]
421 : : end if
422 : :
423 [ + - ]: 876 : if (dtset%useria /= 666) then
424 : 876 : call wfd%copy_cg(ib_v, ik, spin, cg_v)
425 : 876 : call ddkop%apply(ebands%eig(ib_v, ik, spin), npw_k, wfd%nspinor, cg_v, cwaveprj)
426 : : else
427 [ # # ]: 0 : ABI_CHECK(wfd%get_wave_ptr(ib_v, ik, spin, wave_v, msg) == 0, msg)
428 [ # # ]: 0 : ug_v(1:npw_k*nspinor) = wave_v%ug
429 : : end if
430 : :
431 : : ! Loop over bands
432 [ + + ]: 876 : bstop = bmax; if (ds%only_diago) bstop = ib_v
433 [ + + ]: 4996 : do ib_c=ib_v,bstop
434 [ + + ]: 4018 : if (.not. ds%only_diago) then
435 [ + - ]: 3954 : if (distrib_mat(ib_c, ib_v, ik, spin) /= my_rank) cycle
436 : : end if
437 : :
438 [ + - ]: 4894 : if (dtset%useria /= 666) then
439 : 4018 : call wfd%copy_cg(ib_c, ik, spin, cg_c)
440 : 4018 : vv = ddkop%get_braket(ebands%eig(ib_c, ik, spin), istwf_k, npw_k, nspinor, cg_c, mode=ds%mode)
441 : : !if (ib_v == ib_c) vv(2, :) = zero
442 : :
443 [ + + ]: 4018 : if (ds%only_diago) then
444 [ + + ]: 256 : ds%vdiago(:,ib_c,ik,spin) = vv(1, :)
445 : : else
446 [ + + + + ]: 39540 : ds%vmat(:,:,ib_c,ib_v,ik,spin) = vv
447 : : ! Hermitian conjugate
448 [ + + ]: 3954 : if (ib_v /= ib_c) then
449 [ + + ]: 12568 : ds%vmat(1,:,ib_v,ib_c,ik,spin) = vv(1, :)
450 [ + + ]: 12568 : ds%vmat(2,:,ib_v,ib_c,ik,spin) = -vv(2, :)
451 : : end if
452 : : end if
453 : :
454 [ + + ]: 16072 : do idir=1,3
455 [ + + ]: 36162 : ds%dipoles(idir,:,ib_c,ib_v,ik,spin) = vv(:, idir)
456 : : ! Hermitian conjugate
457 [ + + + + ]: 34924 : if (ib_v /= ib_c) ds%dipoles(idir,:,ib_v,ib_c,ik,spin) = [vv(1, idir), -vv(2, idir)]
458 : : end do
459 : :
460 : : else
461 [ # # ]: 0 : ABI_CHECK(wfd%get_wave_ptr(ib_c, ik, spin, wave_c, msg) == 0, msg)
462 [ # # ]: 0 : ug_c(1:npw_k*nspinor) = wave_c%ug
463 : :
464 : : ! Calculate matrix elements of i[H,r] for NC pseudopotentials.
465 [ # # # # : 0 : ihrc = nc_ihr_comm(vkbr, cryst, psps, npw_k, nspinor, istwf_k, dtset%inclvkb, kpt, ug_c, ug_v, kg_k)
# # # # ]
466 : :
467 : : ! HM: 24/07/2018
468 : : ! Transform dipoles to be consistent with results from DFPT
469 : : ! Perturbations with DFPT are along the reciprocal lattice vectors
470 : : ! Perturbations with Commutator are along real space lattice vectors
471 : : ! dot(A, DFPT) = X
472 : : ! dot(B, COMM) = X
473 : : ! B = 2 pi (A^{-1})^T =>
474 : : ! dot(B^T B,COMM) = 2 pi DFPT
475 [ # # # # ]: 0 : vr = (2*pi)*(2*pi)*sum(ihrc(:,:),dim=2)
476 [ # # ]: 0 : vg(1) = dot_product(Cryst%gmet(1,:), vr)
477 [ # # ]: 0 : vg(2) = dot_product(Cryst%gmet(2,:), vr)
478 [ # # ]: 0 : vg(3) = dot_product(Cryst%gmet(3,:), vr)
479 : :
480 : : ! Save matrix elements of i*r in the IBZ
481 [ # # ]: 0 : ds%dipoles(:,1,ib_c,ib_v,ik,spin) = real(vg, kind=dp)
482 [ # # ]: 0 : ds%dipoles(:,1,ib_v,ib_c,ik,spin) = real(vg, kind=dp) ! Hermitian conjugate
483 [ # # ]: 0 : if (ib_v == ib_c) then
484 [ # # ]: 0 : ds%dipoles(:,2,ib_c,ib_v,ik,spin) = zero
485 [ # # ]: 0 : ds%dipoles(:,2,ib_v,ib_c,ik,spin) = zero
486 : : else
487 [ # # ]: 0 : ds%dipoles(:,2,ib_c,ib_v,ik,spin) = aimag(vg)
488 [ # # ]: 0 : ds%dipoles(:,2,ib_v,ib_c,ik,spin) = -aimag(vg) ! Hermitian conjugate
489 : : end if
490 : : end if
491 : :
492 : : end do
493 : : end do
494 : :
495 : : ! Free KB form factors
496 : 102 : call vkbr_free(vkbr)
497 : :
498 [ - + - - : 209 : if (nkpt < 1000 .or. (nkpt > 1000 .and. mod(ik, 200) == 0) .or. ik <= nproc) then
- - - - ]
499 : 102 : write(msg,'(2(a,i0),a)')" k-point [", ik, "/", nkpt, "]"
500 : 102 : call cwtime_report(msg, cpu, wall, gflops)
501 : : end if
502 : :
503 : : end do ! k-points
504 : : end do ! spin
505 : :
506 : 5 : call cwtime_report(msg, cpu_all, wall_all, gflops_all)
507 : :
508 [ - + ]: 5 : ABI_FREE(ug_c)
509 [ - + ]: 5 : ABI_FREE(ug_v)
510 [ - + ]: 5 : ABI_FREE(kg_k)
511 [ - + ]: 5 : ABI_FREE(ihrc)
512 [ + - - + : 5 : ABI_FREE(cwaveprj)
- - - - -
+ ]
513 [ + + ]: 5 : ABI_SFREE(distrib_mat)
514 [ + + ]: 5 : ABI_SFREE(distrib_diago)
515 : :
516 [ + - ]: 5 : if (dtset%useria /= 666) then
517 [ - + ]: 5 : ABI_FREE(cg_c)
518 [ - + ]: 5 : ABI_FREE(cg_v)
519 : 5 : call ddkop%free()
520 : : end if
521 : :
522 : : ! Gather the k-points computed by all processes
523 : 5 : call xmpi_sum_master(ds%dipoles, master, comm, ierr)
524 : :
525 [ + + ]: 5 : if (ds%only_diago) then
526 : 1 : call xmpi_sum_master(ds%vdiago, master, comm, ierr)
527 : : else
528 : 4 : call xmpi_sum_master(ds%vmat, master, comm, ierr)
529 : : end if
530 : :
531 : : ! Write matrix elements to disk.
532 : :
533 : : ! Output EVK file in netcdf format.
534 [ + - + - : 5 : if (my_rank == master .and. write_ncfile .and. dtset%prtevk == 1) then
+ + ]
535 : : ! Have to build hdr on k-grid with info about perturbation.
536 : 2 : call hdr%copy(tmp_hdr)
537 [ + + ]: 8 : tmp_hdr%qptn = zero
538 : :
539 : : !fname = strcat(prefix, "NEW_EVK.nc")
540 : : !call wrtout(ab_out, sjoin("- Writing file: ", fname))
541 : : !NCF_CHECK_MSG(nctk_open_create(ncid, fname, xmpi_comm_self), "Creating EVK.nc file")
542 : : !tmp_hdr%pertcase = 0
543 : : !NCF_CHECK(tmp_hdr%ncwrite(ncid, 43, nc_define=.True.))
544 : : !NCF_CHECK(cryst%ncwrite(ncid))
545 : : !NCF_CHECK(ebands%ncwrite(ncid))
546 : : !if (ds%only_diago) then
547 : : ! ncerr = nctk_def_arrays(ncid, [ &
548 : : ! nctkarr_t('vred_diagonal', "dp", "three, max_number_of_states, number_of_kpoints, number_of_spins")], defmode=.True.)
549 : : !else
550 : : ! ncerr = nctk_def_arrays(ncid, [ nctkarr_t('vred_matrix', "dp", &
551 : : ! "two, three, max_number_of_states, max_number_of_states, number_of_kpoints, number_of_spins")], defmode=.True.)
552 : : !end if
553 : : !NCF_CHECK(ncerr)
554 : : !NCF_CHECK(nctk_set_datamode(ncid))
555 : : !if (ds%only_diago) then
556 : : ! NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "vred_diagonal"), ds%vdiago))
557 : : !else
558 : : ! NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "vred_matrix"), ds%vmat))
559 : : !end if
560 : : !NCF_CHECK(nf90_close(ncid))
561 : :
562 [ + + ]: 8 : do ii=1,3
563 [ + - ]: 6 : fname = strcat(prefix, '_', itoa(ii), "_EVK.nc")
564 : 6 : call wrtout(ab_out, sjoin("- Writing EVK file: ", fname, "for reduced direction:", itoa(ii)))
565 [ - + ]: 6 : NCF_CHECK_MSG(nctk_open_create(ncid, fname, xmpi_comm_self), "Creating EVK.nc file")
566 : 6 : tmp_hdr%pertcase = 3 * cryst%natom + ii
567 [ - + ]: 6 : NCF_CHECK(tmp_hdr%ncwrite(ncid, 43, nc_define=.True.))
568 [ - + ]: 6 : NCF_CHECK(cryst%ncwrite(ncid))
569 [ - + ]: 6 : NCF_CHECK(ebands%ncwrite(ncid))
570 : : ncerr = nctk_def_arrays(ncid, [ &
571 : : nctkarr_t('h1_matrix_elements', "dp", &
572 [ + + ]: 12 : "two, max_number_of_states, max_number_of_states, number_of_kpoints, number_of_spins")], defmode=.True.)
573 [ - + ]: 6 : NCF_CHECK(ncerr)
574 [ - + ]: 6 : NCF_CHECK(nctk_set_datamode(ncid))
575 [ - + ]: 6 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "h1_matrix_elements"), ds%dipoles(ii,:,:,:,:,:)))
576 [ - + ]: 8 : NCF_CHECK(nf90_close(ncid))
577 : : end do
578 : 2 : call tmp_hdr%free()
579 : : end if
580 : :
581 [ + - + - ]: 5 : if (my_rank == master .and. dtset%prtvol > 0) then
582 : 5 : write(ab_out, "(2a)")ch10,"Writing velocity matrix elements (only diagonal terms, real part) for testing purpose:"
583 [ + + ]: 10 : do spin=1,nsppol
584 [ + + ]: 30 : do ik=1,min(nkpt, 4)
585 [ + - ]: 20 : write(ab_out, "(2(a,1x,i0),2x,2a)")"For spin: ", spin, ", ikbz: ", ik, ", kpt: ", trim(ktoa(wfd%kibz(:,ik)))
586 [ + + ]: 192 : do ib_c=bmin,min(bmin+8, bmax)
587 : 192 : write(ab_out, "(3(es16.6,2x))") ds%dipoles(:,1,ib_c,ib_c,ik,spin)
588 : : end do
589 : 25 : write(ab_out,*)""
590 : : !do ib_c=bmin,min(bmin+8, bmax)
591 : : ! write(ab_out, "(a, 6(es16.6,2x))")"Sum_k: ", sum(ds%dipoles(:,:,ib_c,ib_c,:,spin), dim=3) / nkpt
592 : : !end do
593 : : end do
594 : : end do
595 : : end if
596 : :
597 : : ! Free memory
598 : 5 : call wfd%free()
599 : 5 : call ebands%free()
600 : 5 : call cryst%free()
601 : 5 : call hdr%free()
602 : :
603 : : ! Block all procs here so that we know output files are available when code returns.
604 : 5 : call xmpi_barrier(comm)
605 : :
606 [ - + - + : 65 : end subroutine ddkstore_compute_ddk
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - - - -
- - - - -
- - - - -
- + - - -
- - - - -
- - - - -
- - - - -
- - - - -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - - -
- - - - -
- - - - -
- - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- - - - -
- - - - -
- - - - -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + + +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + + +
- + - + +
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - -
- - - - -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + ]
607 : : !!***
608 : :
609 : : !----------------------------------------------------------------------
610 : :
611 : : !!****f* m_ddk/ddkstore_free
612 : : !! NAME
613 : : !! ddkstore_free
614 : : !!
615 : : !! FUNCTION
616 : : !! Free memory
617 : : !!
618 : : !! INPUTS
619 : : !!
620 : : !! SOURCE
621 : :
622 : 5 : subroutine ddkstore_free(self)
623 : :
624 : : !Arguments ------------------------------------
625 : : class(ddkstore_t),intent(inout) :: self
626 : : !************************************************************************
627 : :
628 [ + + ]: 5 : ABI_SFREE(self%vdiago)
629 [ + + ]: 5 : ABI_SFREE(self%vmat)
630 [ + - ]: 5 : ABI_SFREE(self%dipoles)
631 : :
632 : 5 : end subroutine ddkstore_free
633 : : !!***
634 : :
635 : : !----------------------------------------------------------------------
636 : :
637 : : !!****f* m_ddk/ddk_red2car
638 : : !! NAME
639 : : !! ddk_red2car
640 : : !!
641 : : !! FUNCTION
642 : : !! Convert ddk matrix element from reduced coordinates to cartesian coordinates.
643 : : !!
644 : : !! SOURCE
645 : :
646 : 1388 : pure subroutine ddk_red2car(rprimd, vred, vcar)
647 : :
648 : : !Arguments -------------------------------------
649 : : real(dp),intent(in) :: rprimd(3,3)
650 : : real(dp),intent(in) :: vred(2,3)
651 : : real(dp),intent(out) :: vcar(2,3)
652 : :
653 : : !Local variables -------------------------------
654 : : real(dp) :: vtmp(2,3)
655 : : !************************************************************************
656 : :
657 : : ! Go to Cartesian coordinates (same as pmat2cart routine)
658 : : ! V_cart = 1/(2pi) * Rprimd x V_red
659 : : ! where V_red is the derivative computed in the DFPT routines (derivative wrt reduced component).
660 : : vtmp(1,:) = rprimd(:,1)*vred(1,1) &
661 : : +rprimd(:,2)*vred(1,2) &
662 [ + + ]: 5552 : +rprimd(:,3)*vred(1,3)
663 : : vtmp(2,:) = rprimd(:,1)*vred(2,1) &
664 : : +rprimd(:,2)*vred(2,2) &
665 [ + + ]: 5552 : +rprimd(:,3)*vred(2,3)
666 [ + + + + ]: 13880 : vcar = vtmp / two_pi
667 : :
668 : 1388 : end subroutine ddk_red2car
669 : : !!***
670 : :
671 : : !----------------------------------------------------------------------
672 : :
673 : : !!****f* m_ddk/ddkop_init
674 : : !! NAME
675 : : !! ddkop_init
676 : : !!
677 : : !! FUNCTION
678 : : !! Build new object. Use dtset%inclvkb to determine whether the non-local part should be included.
679 : : !!
680 : : !! INPUTS
681 : : !! dtset<dataset_type>=All input variables for this dataset.
682 : : !! cryst<crystal_t>=Crystal structure.
683 : : !! pawtab(ntypat*usepaw)<pawtab_type>=Paw tabulated starting data.
684 : : !! psps<pseudopotential_type>=Variables related to pseudopotentials.
685 : : !! mpi_enreg=information about MPI parallelization
686 : : !! mpw=Maximum number of plane-waves over k-points.
687 : : !! ngfft(18)=contain all needed information about 3D FFT
688 : : !!
689 : : !! SOURCE
690 : :
691 [ + - ]: 76 : subroutine ddkop_init(ddkop, dtset, cryst, pawtab, psps, mpi_enreg, mpw, ngfft)
692 : :
693 : : !Arguments ------------------------------------
694 : : !scalars
695 : : class(ddkop_t),intent(out) :: ddkop
696 : : type(dataset_type),intent(in) :: dtset
697 : : type(crystal_t),intent(in) :: cryst
698 : : type(pseudopotential_type),intent(in) :: psps
699 : : type(MPI_type),target,intent(in) :: mpi_enreg
700 : : integer,intent(in) :: mpw
701 : : !arrays
702 : : integer,intent(in) :: ngfft(18)
703 : : type(pawtab_type),intent(in) :: pawtab(psps%ntypat*psps%usepaw)
704 : :
705 : : !Local variables-------------------------------
706 : : !scalars
707 : : integer,parameter :: cplex1 = 1
708 : : integer :: nfft, mgfft, idir
709 : : ! *************************************************************************
710 : :
711 [ - + ]: 76 : ABI_CHECK(dtset%usepaw == 0, "PAW not tested/implemented!")
712 : :
713 : 76 : ddkop%inclvkb = dtset%inclvkb
714 : 76 : ddkop%usepaw = dtset%usepaw
715 : 76 : ddkop%ipert = cryst%natom + 1
716 : 76 : ddkop%dfpt_sciss = dtset%dfpt_sciss
717 : 76 : ddkop%mpw = mpw
718 : :
719 [ + + + + ]: 988 : ddkop%rprimd = cryst%rprimd
720 : 76 : ddkop%mpi_enreg => mpi_enreg
721 : :
722 : : ! Not used because vlocal1 is not applied.
723 [ + + ]: 304 : nfft = product(ngfft(1:3))
724 [ + + ]: 304 : mgfft = maxval(ngfft(1:3))
725 : :
726 [ + + ]: 304 : do idir=1,3
727 : : ! ==== Initialize most of the Hamiltonian (and derivative) ====
728 : : ! 1) Allocate all arrays and initialize quantities that do not depend on k and spin.
729 : : ! 2) Perform the setup needed for the non-local factors:
730 : : ! * Norm-conserving: Constant kleimann-Bylander energies are copied from psps to gs_hamk.
731 : : ! * PAW: Initialize the overlap coefficients and allocate the Dij coefficients.
732 : : call ddkop%gs_hamkq(idir)%init(psps, pawtab, dtset%nspinor, dtset%nsppol, dtset%nspden, cryst%natom,&
733 : : cryst%typat, cryst%xred, nfft, mgfft, ngfft, cryst%rprimd, dtset%nloalg &
734 : 228 : )
735 : : !paw_ij=paw_ij,comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab,mpi_spintab=mpi_enreg%my_isppoltab,&
736 : : !usecprj=usecprj,ph1d=ph1d,nucdipmom=dtset%nucdipmom,
737 : : ! TODO
738 : : !gpu_option=dtset%gpu_option)
739 : :
740 : : ! Prepare application of the NL part.
741 : 304 : call ddkop%rf_hamkq(idir)%init(cplex1, ddkop%gs_hamkq(idir), ddkop%ipert, has_e1kbsc=.true.)
742 : : end do
743 : :
744 : 76 : end subroutine ddkop_init
745 : : !!***
746 : :
747 : : !----------------------------------------------------------------------
748 : :
749 : : !!****f* m_ddk/ddkop_setup_spin_kpoint
750 : : !! NAME
751 : : !! ddkop_setup_spin_kpoint
752 : : !!
753 : : !! FUNCTION
754 : : !! Prepare internal tables that depend on k-point/spin
755 : : !!
756 : : !! INPUTS
757 : : !! dtset<dataset_type>=All input variables for this dataset.
758 : : !! cryst<crystal_t>=Crystal structure.
759 : : !! psps<pseudopotential_type>=Variables related to pseudopotentials.
760 : : !! spin: spin index
761 : : !! kpoint(3): K-point in reduced coordinates.
762 : : !! istwkf_k: defines storage of wavefunctions for this k-point
763 : : !! npw_k: Number of planewaves.
764 : : !! kg_k(3,npw_k)=reduced planewave coordinates.
765 : : !!
766 : : !! SOURCE
767 : :
768 : 433 : subroutine ddkop_setup_spin_kpoint(ddkop, dtset, cryst, psps, spin, kpoint, istwf_k, npw_k, kg_k)
769 : :
770 : : !Arguments ------------------------------------
771 : : !scalars
772 : : class(ddkop_t),intent(inout) :: ddkop
773 : : type(dataset_type),intent(in) :: dtset
774 : : type(crystal_t),intent(in) :: cryst
775 : : type(pseudopotential_type),intent(in) :: psps
776 : : integer,intent(in) :: spin, npw_k, istwf_k
777 : : !arrays
778 : : integer,intent(in) :: kg_k(3,npw_k)
779 : : real(dp),intent(in) :: kpoint(3)
780 : :
781 : : !Local variables-------------------------------
782 : : !scalars
783 : : integer,parameter :: nkpt1=1, nsppol1=1
784 : 433 : type(mpi_type) :: mpienreg_seq
785 : : !arrays
786 : : integer :: npwarr(nkpt1), dummy_nband(nkpt1*nsppol1)
787 : : integer :: idir, nkpg, nkpg1, useylmgr1, optder !, nylmgr1
788 : 433 : real(dp),allocatable :: ylm_k(:,:),ylmgr1_k(:,:,:)
789 : : !************************************************************************
790 : :
791 [ - + ]: 433 : ABI_CHECK_ILEQ(npw_k, ddkop%mpw, "npw_k > mpw!")
792 [ + + ]: 1732 : ddkop%kpoint = kpoint
793 : :
794 : : ! Set up the spherical harmonics (Ylm) at k+q if useylm = 1
795 : 433 : useylmgr1 = 0; optder = 0
796 [ + + ]: 433 : if (psps%useylm == 1) then
797 : 18 : useylmgr1 = 1; optder = 1
798 : : end if
799 : :
800 [ + + - + : 2213314 : ABI_RECALLOC(ddkop%gh1c, (2, npw_k*dtset%nspinor, 3))
+ - + - +
- + - - +
- + - + +
+ + + +
+ ]
801 [ + + - + : 2212448 : ABI_RECALLOC(ddkop%gs1c, (2, npw_k*dtset%nspinor, 3))
+ - - + -
+ - + + +
+ + + + ]
802 : :
803 [ - + + - : 1732 : ABI_MALLOC(ylm_k, (npw_k, psps%mpsang**2 * psps%useylm))
+ + + - +
- + + - +
- + ]
804 [ - + + - : 2165 : ABI_MALLOC(ylmgr1_k, (npw_k, 3+6*(optder/2), psps%mpsang**2*psps%useylm*useylmgr1))
+ - - + +
- + + + -
+ - + + -
+ - + ]
805 : :
806 [ + + ]: 433 : if (psps%useylm == 1) then
807 : : ! Fake MPI_type for sequential part. dummy_nband and nsppol1 are not used in sequential mode.
808 : 18 : call initmpi_seq(mpienreg_seq)
809 [ + + ]: 36 : dummy_nband = 0; npwarr = npw_k
810 : : call initylmg(cryst%gprimd, kg_k, kpoint, nkpt1, mpienreg_seq, psps%mpsang, npw_k, dummy_nband, nkpt1, &
811 : 18 : npwarr, nsppol1, optder, cryst%rprimd, ylm_k, ylmgr1_k)
812 : 18 : call destroy_mpi_enreg(mpienreg_seq)
813 : : end if
814 : :
815 [ + + ]: 1732 : do idir=1,3
816 : 1299 : call ddkop%htg(idir)%free()
817 : :
818 : : ! Continue to initialize the Hamiltonian
819 : 1299 : call ddkop%gs_hamkq(idir)%load_spin(spin, with_nonlocal=.true.)
820 : 1299 : call ddkop%rf_hamkq(idir)%load_spin(spin, with_nonlocal=.true.)
821 : :
822 : : ! We need ffnl1 and dkinpw for 3 dirs. Note that the Hamiltonian objects use pointers to keep a reference
823 : : ! to the output results of this routine.
824 : : ! This is the reason why we need to store the targets in ddkop%htg
825 : : call getgh1c_setup(ddkop%gs_hamkq(idir), ddkop%rf_hamkq(idir), dtset, psps, kpoint, kpoint, idir, ddkop%ipert, & ! In
826 : : cryst%natom, cryst%rmet, cryst%gprimd, cryst%gmet, istwf_k, npw_k, npw_k, & ! In
827 : : useylmgr1, kg_k, ylm_k, kg_k, ylm_k, ylmgr1_k, & ! In
828 : : ddkop%htg(idir)%dkinpw, nkpg, nkpg1, ddkop%htg(idir)%kpg_k, ddkop%htg(idir)%kpg1_k, & ! Out
829 : : ddkop%htg(idir)%kinpw1, ddkop%htg(idir)%ffnlk, ddkop%htg(idir)%ffnl1, & ! Out
830 [ - + - + : 1732 : ddkop%htg(idir)%ph3d, ddkop%htg(idir)%ph3d1) ! Out
- + - + ]
831 : : end do
832 : :
833 [ - + ]: 433 : ABI_FREE(ylm_k)
834 [ - + ]: 433 : ABI_FREE(ylmgr1_k)
835 : :
836 [ - + - + : 433 : end subroutine ddkop_setup_spin_kpoint
- + - + -
+ - + - +
- + - + -
+ - + ]
837 : : !!***
838 : :
839 : : !----------------------------------------------------------------------
840 : :
841 : : !!****f* m_ddk/ddkop_apply
842 : : !! NAME
843 : : !! ddkop_apply
844 : : !!
845 : : !! FUNCTION
846 : : !! Apply velocity operator dH/dk to wavefunction in G-space. Store results in object.
847 : : !!
848 : : !! INPUTS
849 : : !! eig0nk: Eigenvalue associated to the wavefunction.
850 : : !! npw_k: Number of planewaves.
851 : : !! nspinor: Number of spinor components.
852 : : !! cwave(2,npw_k*nspinor)=input wavefunction in reciprocal space
853 : : !! cwaveprj(natom,nspinor*usecprj)=<p_lmn|C> coefficients for wavefunction |C> (and 1st derivatives)
854 : : !! if not allocated or size=0, they are locally computed (and not sorted)!!
855 : : !!
856 : : !! SIDE EFFECTS
857 : : !! Stores:
858 : : !! gh1c(2,npw_k*nspinor)= <G|H^(1)|C> or <G|H^(1)-lambda.S^(1)|C> on the k+q sphere
859 : : !! (only kinetic+non-local parts if optlocal=0)
860 : : !!
861 : : !! SOURCE
862 : :
863 [ - + ]: 2264 : subroutine ddkop_apply(ddkop, eig0nk, npw_k, nspinor, cwave, cwaveprj)
864 : :
865 : : !Arguments ------------------------------------
866 : : !scalars
867 : : class(ddkop_t),target,intent(inout) :: ddkop
868 : : integer,intent(in) :: npw_k, nspinor
869 : : real(dp),intent(in) :: eig0nk
870 : : !arrays
871 : : real(dp),intent(inout) :: cwave(2,npw_k*nspinor)
872 : : type(pawcprj_type),intent(inout) :: cwaveprj(:,:)
873 : :
874 : : !Local variables-------------------------------
875 : : !scalars
876 : : integer,parameter :: berryopt0 = 0, optlocal0 = 0, tim_getgh1c = 1, usevnl0 = 0, opt_gvnlx1 = 0
877 : : integer :: idir, sij_opt, ispinor, ipws, ipw, optnl
878 : : real(dp) :: eshift
879 : : !arrays
880 : : real(dp) :: grad_berry(2,(berryopt0/4)), gvnlx1(2,usevnl0)
881 : : !************************************************************************
882 : :
883 : 2264 : ddkop%eig0nk = eig0nk
884 : :
885 [ + + ]: 2264 : if (ddkop%inclvkb /= 0) then
886 : : !if (.True.) then
887 : : ! optlocal0 = 0: local part of H^(1) is not computed in gh1c=<G|H^(1)|C>
888 : : ! optnl = 2: non-local part of H^(1) is totally computed in gh1c=<G|H^(1)|C>
889 : : ! opt_gvnlx1 = option controlling the use of gvnlx1 array:
890 : 2164 : optnl = 2; if (ddkop%inclvkb == 0) optnl = 0
891 : :
892 : 2164 : eshift = ddkop%eig0nk - ddkop%dfpt_sciss
893 [ + + ]: 8656 : do idir=1,3
894 : 6492 : sij_opt = ddkop%gs_hamkq(idir)%usepaw
895 : : call getgh1c(berryopt0, cwave, cwaveprj, ddkop%gh1c(:,:,idir), &
896 : : grad_berry, ddkop%gs1c(:,:,idir), ddkop%gs_hamkq(idir), gvnlx1, idir, ddkop%ipert, [eshift], ddkop%mpi_enreg, 1, optlocal0, &
897 [ + + ]: 15148 : optnl, opt_gvnlx1, ddkop%rf_hamkq(idir), sij_opt, tim_getgh1c, usevnl0)
898 : : end do
899 : :
900 : : else
901 : : ! FIXME: optnl 0 with DDK does not work as expected.
902 : : ! So I treat the kinetic term explicitly without calling getgh1c.
903 [ + + ]: 400 : do idir=1,3
904 : 100 : associate (kinpw1 => ddkop%gs_hamkq(idir)%kinpw_kp, dkinpw => ddkop%rf_hamkq(idir)%dkinpw_k)
905 [ + + ]: 600 : do ispinor=1,nspinor
906 [ + + ]: 346380 : do ipw=1,npw_k
907 : 345780 : ipws = ipw + npw_k*(ispinor-1)
908 [ + - ]: 346080 : if (kinpw1(ipw) < huge(zero)*1.d-11) then
909 [ + + ]: 1037340 : ddkop%gh1c(:,ipws,idir) = dkinpw(ipw) * cwave(:,ipws)
910 : : else
911 [ # # ]: 0 : ddkop%gh1c(:,ipws,idir) = zero
912 : : end if
913 : : end do
914 : : end do
915 : : end associate
916 : : end do ! idir
917 : : end if
918 : :
919 : 2264 : end subroutine ddkop_apply
920 : : !!***
921 : :
922 : : !----------------------------------------------------------------------
923 : :
924 : : !!****f* m_ddk/ddkop_get_braket
925 : : !! NAME
926 : : !! ddkop_get_braket
927 : : !!
928 : : !! FUNCTION
929 : : !! Compute diagonal matrix element in Cartesian coordinates.
930 : : !!
931 : : !! INPUTS
932 : : !! eig0mk: Eigenvalue associated to the "bra" wavefunction
933 : : !! istwkf_k: defines storage of wavefunctions for this k-point
934 : : !! npw_k: Number of planewaves.
935 : : !! nspinor: Number of spinor components.
936 : : !! brag(2,npw_k*nspinor)=input wavefunction in reciprocal space
937 : : !!
938 : : !! SOURCE
939 : :
940 : 10812 : function ddkop_get_braket(ddkop, eig0mk, istwf_k, npw_k, nspinor, brag, mode) result(vk)
941 : :
942 : : !Arguments ------------------------------------
943 : : !scalars
944 : : class(ddkop_t),intent(in) :: ddkop
945 : : integer,intent(in) :: istwf_k, npw_k, nspinor
946 : : real(dp),intent(in) :: eig0mk
947 : : character(len=*),optional,intent(in) :: mode
948 : : !arrays
949 : : real(dp),intent(in) :: brag(2*npw_k*nspinor)
950 : : real(dp) :: vk(2,3)
951 : :
952 : : !Local variables-------------------------------
953 : : !scalars
954 : : integer :: idir
955 : : real(dp) :: doti
956 : : !arrays
957 : : real(dp) :: dotarr(2), vk_red(2, 3)
958 : : character(len=50) :: my_mode
959 : : !************************************************************************
960 : :
961 [ + - ]: 5406 : if (ddkop%usepaw == 0) then
962 : : ! <u_(iband,k+q)^(0)|H_(k+q,k)^(1)|u_(jband,k)^(0)> (NC psps)
963 [ + + ]: 21624 : do idir=1,3
964 : 16218 : dotarr = cg_zdotc(npw_k * nspinor, brag, ddkop%gh1c(:,:,idir))
965 [ + + ]: 16218 : if (istwf_k > 1) then
966 : : !dum = two * j_dpc * AIMAG(dum); if (vkbr%istwfk==2) dum = dum - j_dpc * AIMAG(gamma_term)
967 : 72 : doti = two * dotarr(2)
968 [ + + + - ]: 72 : if (istwf_k == 2 .and. ddkop%mpi_enreg%me_g0 == 1) then
969 : : ! nspinor always 1
970 : : ! TODO: Recheck this part but it should be ok.
971 : 24 : doti = doti - (brag(1) * ddkop%gh1c(2,1,idir) - brag(2) * ddkop%gh1c(1,1,idir))
972 : : end if
973 : 72 : dotarr(2) = doti; dotarr(1) = zero
974 : : end if
975 [ + + ]: 54060 : vk(:, idir) = dotarr
976 : : end do
977 : : else
978 : 0 : ABI_ERROR("PAW Not Implemented")
979 : : ! <u_(iband,k+q)^(0)|H_(k+q,k)^(1)-(eig0_k+eig0_k+q)/2.S^(1)|u_(jband,k)^(0)> (PAW)
980 : : ! eshiftkq = half * (eig0mk - ddkop%eig0nk)
981 : : ABI_UNUSED(eig0mk)
982 : : end if
983 : :
984 [ + - - + ]: 5406 : my_mode = "cart"; if (present(mode)) my_mode = mode
985 : : select case (mode)
986 : : case ("cart")
987 [ + + + + ]: 13880 : vk_red = vk
988 : 1388 : call ddk_red2car(ddkop%rprimd, vk_red, vk)
989 : : case ("reduced")
990 : 0 : continue
991 : : case default
992 [ + - + ]: 5406 : ABI_ERROR(sjoin("Invalid vaue for mode:", mode))
993 : : end select
994 : :
995 [ - + ]: 5406 : end function ddkop_get_braket
996 : : !!***
997 : :
998 : : !----------------------------------------------------------------------
999 : :
1000 : : !!****f* m_ddk/ddkop_get_ihr_comm
1001 : : !! NAME
1002 : : !! ddkop_get_ihr_comm
1003 : : !!
1004 : : !! FUNCTION
1005 : : !! Compute matrix element in Cartesian coordinates.
1006 : : !!
1007 : : !! INPUTS
1008 : : !! eig0mk: Eigenvalue associated to the "bra" wavefunction
1009 : : !! istwkf_k: defines storage of wavefunctions for this k-point
1010 : : !! npw_k: Number of planewaves.
1011 : : !! nspinor: Number of spinor components.
1012 : : !! brag(2,npw_k*nspinor)=input wavefunction in reciprocal space
1013 : : !!
1014 : : !! SOURCE
1015 : :
1016 : 0 : subroutine ddkop_get_ihr_comm(ddkop, cryst, eig0mk, istwf_k, npw_k, nspinor, brag, new_rhotwx)
1017 : :
1018 : : !Arguments ------------------------------------
1019 : : !scalars
1020 : : class(ddkop_t),intent(in) :: ddkop
1021 : : type(crystal_t),intent(in) :: cryst
1022 : : integer,intent(in) :: istwf_k, npw_k, nspinor
1023 : : real(dp),intent(in) :: eig0mk
1024 : : !arrays
1025 : : real(dp),intent(in) :: brag(2,npw_k*nspinor)
1026 : : complex(gwp),intent(out) :: new_rhotwx(3, nspinor**2)
1027 : :
1028 : : !Local variables-------------------------------
1029 : : !scalars
1030 : : integer :: idir, iab
1031 : : real(dp) :: doti
1032 : : !arrays
1033 : : integer :: spinor_pad(2,4), spad1, spad2 !, rtmp(2)
1034 : : real(dp) :: dotarr(2)
1035 : : ! real(dp), allocatable :: cg2_dp(:,:), ddk_ug1(:,:,:)
1036 : : !************************************************************************
1037 : :
1038 [ # # ]: 0 : spinor_pad = reshape([0, 0, npw_k, npw_k, 0, npw_k, npw_k, 0], [2, 4])
1039 : :
1040 [ # # ]: 0 : if (ddkop%usepaw == 0) then
1041 : : ! <u_(iband,k+q)^(0)|H_(k+q,k)^(1)|u_(jband,k)^(0)> (NC psps)
1042 [ # # ]: 0 : do iab=1,nspinor**2
1043 : 0 : spad1 = spinor_pad(1,iab); spad2 = spinor_pad(2,iab)
1044 [ # # ]: 0 : do idir=1,3
1045 : 0 : dotarr = cg_zdotc(npw_k, brag(:,spad2+1), ddkop%gh1c(:,spad1+1,idir))
1046 [ # # ]: 0 : if (istwf_k > 1) then
1047 : 0 : doti = two * dotarr(2)
1048 [ # # # # ]: 0 : if (istwf_k == 2 .and. ddkop%mpi_enreg%me_g0 == 1) then
1049 : : ! nspinor always 1
1050 : : ! TODO: Recheck this part but it should be ok.
1051 : 0 : doti = doti - (brag(1,1) * ddkop%gh1c(2,1,idir) - brag(2,1) * ddkop%gh1c(1,1,idir))
1052 : : end if
1053 : 0 : dotarr(2) = doti; dotarr(1) = zero
1054 : : end if
1055 : 0 : new_rhotwx(idir, iab) = dotarr(1) + j_dpc * dotarr(2)
1056 : : end do
1057 : : end do ! iab
1058 : : else
1059 : 0 : ABI_ERROR("PAW Not Implemented")
1060 : : ! <u_(iband,k+q)^(0)|H_(k+q,k)^(1)-(eig0_k+eig0_k+q)/2.S^(1)|u_(jband,k)^(0)> (PAW)
1061 : : ! eshiftkq = half * (eig0mk - ddkop%eig0nk)
1062 : : ABI_UNUSED(eig0mk)
1063 : : end if
1064 : :
1065 : : ! HM: 24/07/2018
1066 : : ! Transform dipoles to be consistent with results from DFPT
1067 : : ! Perturbations with DFPT are along the reciprocal lattice vectors
1068 : : ! Perturbations with commutator are along real space lattice vectors
1069 : : ! dot(A, DFPT) = X
1070 : : ! dot(B, COMM) = X
1071 : : ! B = 2 pi (A^{-1})^T => dot(B^T B,COMM) = 2 pi DFPT
1072 : :
1073 [ # # ]: 0 : do iab=1,nspinor**2
1074 [ # # # # : 0 : new_rhotwx(:, iab) = matmul(cryst%rmet, new_rhotwx(:, iab)) / (two_pi ** 2)
# # # # #
# ]
1075 : : end do
1076 : :
1077 : 0 : end subroutine ddkop_get_ihr_comm
1078 : : !!***
1079 : :
1080 : : !----------------------------------------------------------------------
1081 : :
1082 : : !!****f* m_ddk/ddkop_get_vdiag
1083 : : !! NAME
1084 : : !! ddkop_get_vdiag
1085 : : !!
1086 : : !! FUNCTION
1087 : : !! Simplified interface to compute the diagonal matrix element of the velocity operator in cartesian coords.
1088 : : !!
1089 : : !! INPUTS
1090 : : !!
1091 : : !! SOURCE
1092 : :
1093 [ - + ]: 1388 : function ddkop_get_vdiag(ddkop, eig0nk, istwf_k, npw_k, nspinor, cwave, cwaveprj, mode) result(vk)
1094 : :
1095 : : !Arguments ------------------------------------
1096 : : !scalars
1097 : : class(ddkop_t),intent(inout) :: ddkop
1098 : : integer,intent(in) :: istwf_k, npw_k, nspinor
1099 : : real(dp),intent(in) :: eig0nk
1100 : : character(len=*),optional,intent(in) :: mode
1101 : : !arrays
1102 : : real(dp),intent(inout) :: cwave(2,npw_k*nspinor)
1103 : : type(pawcprj_type),intent(inout) :: cwaveprj(:,:)
1104 : : real(dp) :: vk(3)
1105 : :
1106 : : !Local variables-------------------------------
1107 : : character(len=50) :: my_mode
1108 : : !arrays
1109 : : real(dp) :: cvk(2, 3)
1110 : : !************************************************************************
1111 : :
1112 [ - + - - ]: 1388 : my_mode = "cart"; if (present(mode)) my_mode = mode
1113 : 1388 : call ddkop%apply(eig0nk, npw_k, nspinor, cwave, cwaveprj)
1114 [ + + + + ]: 13880 : cvk = ddkop%get_braket(eig0nk, istwf_k, npw_k, nspinor, cwave, mode=my_mode)
1115 [ + + ]: 5552 : vk = cvk(1, :)
1116 : :
1117 [ - + ]: 1388 : end function ddkop_get_vdiag
1118 : : !!***
1119 : :
1120 : : !!****f* m_ddk/ddkop_get_vnondiag
1121 : : !! NAME
1122 : : !! ddkop_get_vdiag
1123 : : !!
1124 : : !! FUNCTION
1125 : : !! Simplified interface to compute the off-diagonal matrix elemente of the velocity operator in cartesian coords.
1126 : : !!
1127 : : !! INPUTS
1128 : : !!
1129 : : !! SOURCE
1130 : :
1131 [ # # ]: 0 : function ddkop_get_vnondiag(ddkop, eig0nk_bra, istwf_k, npw_k, nspinor, cwave_bra, cwave_ket, cwaveprj, mode) result(cvk)
1132 : :
1133 : : !Arguments ------------------------------------
1134 : : !scalars
1135 : : class(ddkop_t),intent(inout) :: ddkop
1136 : : integer,intent(in) :: istwf_k, npw_k, nspinor
1137 : : real(dp),intent(in) :: eig0nk_bra
1138 : : character(len=*),optional,intent(in) :: mode
1139 : : !arrays
1140 : : real(dp),intent(inout) :: cwave_bra(2,npw_k*nspinor),cwave_ket(2,npw_k*nspinor)
1141 : : type(pawcprj_type),intent(inout) :: cwaveprj(:,:)
1142 : : real(dp) :: cvk(2,3)
1143 : :
1144 : : !Local variables-------------------------------
1145 : : character(len=50) :: my_mode
1146 : : !************************************************************************
1147 : :
1148 [ # # # # ]: 0 : my_mode = "cart"; if (present(mode)) my_mode = mode
1149 : 0 : call ddkop%apply(eig0nk_bra, npw_k, nspinor, cwave_ket, cwaveprj)
1150 [ # # # # ]: 0 : cvk = ddkop%get_braket(eig0nk_bra, istwf_k, npw_k, nspinor, cwave_bra, mode=my_mode)
1151 : :
1152 [ # # ]: 0 : end function ddkop_get_vnondiag
1153 : : !!***
1154 : :
1155 : : !----------------------------------------------------------------------
1156 : :
1157 : : !!****f* m_ddk/ddkop_free
1158 : : !! NAME
1159 : : !! ddkop_free
1160 : : !!
1161 : : !! FUNCTION
1162 : : !! Free memory
1163 : : !!
1164 : : !! SOURCE
1165 : :
1166 : 76 : subroutine ddkop_free(ddkop)
1167 : :
1168 : : !Arguments ------------------------------------
1169 : : !scalars
1170 : : class(ddkop_t),intent(inout) :: ddkop
1171 : :
1172 : : !Local variables-------------------------------
1173 : : integer :: idir
1174 : : !************************************************************************
1175 : :
1176 [ + - ]: 76 : ABI_SFREE(ddkop%gh1c)
1177 [ + - ]: 76 : ABI_SFREE(ddkop%gs1c)
1178 : :
1179 [ + + ]: 304 : do idir=1,3
1180 : 228 : call ddkop%gs_hamkq(idir)%free()
1181 : 228 : call ddkop%htg(idir)%free()
1182 : 304 : call ddkop%rf_hamkq(idir)%free()
1183 : : end do
1184 : :
1185 : 76 : ddkop%mpi_enreg => null()
1186 : :
1187 : 76 : end subroutine ddkop_free
1188 : : !!***
1189 : :
1190 : : !----------------------------------------------------------------------
1191 : :
1192 : : !!****f* m_ddk/ham_targets_free
1193 : : !! NAME
1194 : : !!
1195 : : !! FUNCTION
1196 : : !!
1197 : : !! INPUTS
1198 : : !!
1199 : : !! SOURCE
1200 : :
1201 : 1527 : subroutine ham_targets_free(ddkop)
1202 : :
1203 : : !Arguments ------------------------------------
1204 : : class(ham_targets_t),intent(inout) :: ddkop
1205 : : !************************************************************************
1206 : :
1207 [ + + ]: 1527 : ABI_SFREE(ddkop%ffnlk)
1208 [ + + ]: 1527 : ABI_SFREE(ddkop%ffnl1)
1209 [ + + ]: 1527 : ABI_SFREE(ddkop%kpg_k)
1210 [ + + ]: 1527 : ABI_SFREE(ddkop%kpg1_k)
1211 [ + + ]: 1527 : ABI_SFREE(ddkop%dkinpw)
1212 [ + + ]: 1527 : ABI_SFREE(ddkop%kinpw1)
1213 [ + + ]: 1527 : ABI_SFREE(ddkop%ph3d)
1214 [ - + ]: 1527 : ABI_SFREE(ddkop%ph3d1)
1215 : :
1216 : 1527 : end subroutine ham_targets_free
1217 : : !!***
1218 : :
1219 : : !----------------------------------------------------------------------
1220 : :
1221 [ - + - - : 5624 : end module m_ddk
+ + - + +
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + +
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + +
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + +
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + +
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + +
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + +
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + +
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + +
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + +
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + +
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + +
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + +
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + +
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + +
+ - + - +
+ + - + -
+ + + - +
- + - + -
+ - + - +
- + - + +
+ - + - +
- + - + -
+ - + - +
- + + + -
+ - + - +
- + - + -
+ - + - +
+ + - + -
+ - + - +
- + - + -
+ - + + +
- + - + -
+ - + - +
- + - + -
+ + + - +
- + - + -
+ - + - +
- + - + +
+ - + - +
- + - + -
+ - + - +
- + + + -
+ - + - +
- + - + -
+ - + - +
+ - - + -
+ ]
1222 : : !!***
|
