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 3 : 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 3 : type(wfd_t) :: wfd
233 3 : type(vkbr_t) :: vkbr
234 3 : type(ebands_t) :: ebands
235 156 : type(crystal_t) :: cryst
236 3 : type(hdr_type) :: tmp_hdr, hdr
237 30 : type(ddkop_t) :: ddkop
238 : type(wave_t),pointer :: wave_v, wave_c
239 : !arrays
240 3 : integer,allocatable :: distrib_mat(:,:,:,:), distrib_diago(:,:,:), nband(:,:), kg_k(:,:)
241 6 : logical,allocatable :: bks_mask(:,:,:), keep_ur(:,:,:)
242 : real(dp) :: kpt(3), vv(2, 3)
243 3 : real(dp),allocatable :: cg_c(:,:), cg_v(:,:)
244 : complex(dp) :: vg(3), vr(3)
245 3 : complex(gwp),allocatable :: ihrc(:,:), ug_c(:), ug_v(:)
246 3 : type(pawcprj_type),allocatable :: cwaveprj(:,:)
247 : !************************************************************************
248 :
249 3 : my_rank = xmpi_comm_rank(comm); nproc = xmpi_comm_size(comm)
250 :
251 9 : if (my_rank == master) call wrtout([std_out, ab_out], " Computation of velocity matrix elements (DDK)", newlines=1)
252 :
253 3 : ABI_CHECK(psps%usepaw == 0, "PAW not implemented")
254 :
255 : ! Get ebands and hdr from WFK file.
256 3 : ebands = wfk_read_ebands(wfk_path, comm, out_hdr=hdr)
257 3 : cryst = hdr%get_crystal()
258 :
259 : ! Extract important dimensions from hdr%
260 3 : nkpt = hdr%nkpt
261 3 : nsppol = hdr%nsppol
262 3 : nspinor = hdr%nspinor
263 3 : 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 3 : if (ds%bmax == -1) ds%bmax = mband
269 :
270 3 : 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 3 : bmin = ds%bmin; bmax = ds%bmax
275 3 : nbcalc = bmax - bmin + 1
276 3 : write_ncfile = len_trim(prefix) > 0
277 3 : 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 3 : if (my_rank == master) then
283 3 : write(ab_out, "(a)")" Parameters extracted from the Abinit header:"
284 3 : write(ab_out, "(a, f5.1)") ' ecut: ', hdr%ecut
285 3 : write(ab_out, "(a, i0)") ' nkpt: ', nkpt
286 3 : write(ab_out, "(a, i0)") ' mband: ', mband
287 3 : write(ab_out, "(a, i0)") ' nsppol: ', nsppol
288 3 : write(ab_out, "(a, i0)") ' nspinor: ', nspinor
289 3 : write(ab_out, "(a, i0)") ' useylm: ', dtset%useylm
290 3 : write(ab_out, "(a, i0)") ' inclvkb: ', dtset%inclvkb
291 3 : write(ab_out, "(2(a, i0))")' bmin: ', bmin, ", bmax: ", bmax
292 3 : if (ds%only_diago) then
293 0 : write(ab_out, "(a)")' Computing diagonal matrix elements only'
294 : else
295 3 : write(ab_out, "(a)")' Computing diagonal and off-diagonal matrix elements'
296 : end if
297 3 : write(ab_out, "(2(a, i0))")' Between band index bmin: ', bmin, ", bmax: ", bmax
298 3 : write(ab_out, "(a)")""
299 : end if
300 :
301 : ! Create distribution of the wavefunctions mask.
302 12 : ABI_MALLOC(nband, (nkpt, nsppol))
303 15 : ABI_MALLOC(keep_ur, (mband, nkpt, nsppol))
304 12 : ABI_MALLOC(bks_mask, (mband, nkpt, nsppol))
305 702 : keep_ur = .false.; bks_mask = .false.; nband = mband
306 :
307 3 : if (ds%only_diago) then
308 : ! Distribute k-points, spin and (b, b) diagonal over MPI processors.
309 0 : ABI_MALLOC(distrib_diago, (bmin:bmax, nkpt, nsppol))
310 0 : distrib_diago = -1
311 :
312 : ! Create bks_mask to load the wavefunctions.
313 0 : ii = 0
314 0 : do spin=1,nsppol
315 0 : do ik=1,nkpt
316 0 : do ib_v=bmin,bmax
317 0 : ii = ii + 1; if (mod(ii, nproc) /= my_rank) cycle ! MPI parallelism.
318 0 : distrib_diago(ib_v, ik, spin) = my_rank
319 0 : bks_mask(ib_v, ik, spin) = .true.
320 : end do
321 : end do
322 : end do
323 0 : 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 18 : ABI_MALLOC(distrib_mat, (bmin:bmax, bmin:bmax, nkpt, nsppol))
328 3 : call xmpi_distab(nproc, distrib_mat)
329 :
330 : ! Create bks_mask to load the wavefunctions
331 6 : do spin=1,nsppol
332 36 : do ik=1,nkpt
333 : ! Loop over v bands
334 333 : do ib_v=bmin,bmax
335 : ! Loop over c bands
336 3330 : do ib_c=bmin,bmax
337 3300 : if (distrib_mat(ib_c, ib_v, ik, spin) == my_rank) then
338 3000 : bks_mask(ib_v, ik, spin) = .true.
339 3000 : bks_mask(ib_c, ik, spin) = .true.
340 : end if
341 : end do
342 : end do
343 : end do
344 : end do
345 :
346 3336 : 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 3 : ngfftc, dtset%nloalg, dtset%prtvol, dtset%pawprtvol, comm)
353 :
354 3 : ABI_FREE(bks_mask)
355 3 : ABI_FREE(keep_ur)
356 3 : ABI_FREE(nband)
357 :
358 6 : call wfd%print([std_out], header="Wavefunctions on the k-points grid")
359 :
360 : ! Read wavefunctions from WFK file.
361 3 : call wfd%read_wfk(wfk_path, iomode_from_fname(wfk_path))
362 :
363 : ! Allocate workspace arrays
364 33 : mpw = maxval(wfd%npwarr)
365 9 : ABI_MALLOC(kg_k, (3, mpw))
366 9 : ABI_MALLOC(ug_c, (mpw*nspinor))
367 6 : ABI_MALLOC(ug_v, (mpw*nspinor))
368 3 : if (dtset%useria /= 666) then
369 9 : ABI_MALLOC(cg_c, (2, mpw*nspinor))
370 6 : ABI_MALLOC(cg_v, (2, mpw*nspinor))
371 : end if
372 :
373 3 : ABI_MALLOC(cwaveprj, (0, 0))
374 27351 : ABI_CALLOC(ds%dipoles, (3, 2, bmin:bmax, bmin:bmax, nkpt, nsppol))
375 9 : ABI_MALLOC(ihrc, (3, nspinor**2))
376 :
377 3 : if (ds%only_diago) then
378 0 : ABI_CALLOC(ds%vdiago, (3, bmin:bmax, nkpt, nsppol))
379 : else
380 30348 : ABI_CALLOC(ds%vmat, (2, 3, bmin:bmax, bmin:bmax, nkpt, nsppol))
381 : end if
382 :
383 3 : if (dtset%useria /= 666) then
384 3 : 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 3 : call cwtime(cpu_all, wall_all, gflops_all, "start")
389 :
390 6 : do spin=1,nsppol
391 36 : do ik=1,nkpt
392 :
393 : ! Only do a subset a k-points
394 30 : if (ds%only_diago) then
395 0 : if (all(distrib_diago(:, ik, spin) /= my_rank)) cycle
396 : else
397 30 : if (all(distrib_mat(bmin:bmax, bmin:bmax, ik, spin) /= my_rank)) cycle
398 : end if
399 30 : call cwtime(cpu, wall, gflops, "start")
400 :
401 30 : nband_k = wfd%nband(ik, spin)
402 30 : istwf_k = wfd%istwfk(ik)
403 30 : npw_k = wfd%npwarr(ik)
404 120 : kpt = wfd%kibz(:,ik)
405 138342 : kg_k(:,1:npw_k) = wfd%kdata(ik)%kg_k
406 :
407 30 : if (dtset%useria /= 666) then
408 30 : 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 330 : do ib_v=bmin,bmax
417 300 : if (ds%only_diago) then
418 0 : if (distrib_diago(ib_v,ik,spin) /= my_rank) cycle
419 : else
420 300 : if (all(distrib_mat(:,ib_v,ik,spin) /= my_rank)) cycle
421 : end if
422 :
423 300 : if (dtset%useria /= 666) then
424 300 : call wfd%copy_cg(ib_v, ik, spin, cg_v)
425 300 : 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 300 : bstop = bmax; if (ds%only_diago) bstop = ib_v
433 1980 : do ib_c=ib_v,bstop
434 1650 : if (.not. ds%only_diago) then
435 1650 : if (distrib_mat(ib_c, ib_v, ik, spin) /= my_rank) cycle
436 : end if
437 :
438 1950 : if (dtset%useria /= 666) then
439 1650 : call wfd%copy_cg(ib_c, ik, spin, cg_c)
440 1650 : 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 1650 : if (ds%only_diago) then
444 0 : ds%vdiago(:,ib_c,ik,spin) = vv(1, :)
445 : else
446 16500 : ds%vmat(:,:,ib_c,ib_v,ik,spin) = vv
447 : ! Hermitian conjugate
448 1650 : if (ib_v /= ib_c) then
449 5400 : ds%vmat(1,:,ib_v,ib_c,ik,spin) = vv(1, :)
450 5400 : ds%vmat(2,:,ib_v,ib_c,ik,spin) = -vv(2, :)
451 : end if
452 : end if
453 :
454 6600 : do idir=1,3
455 14850 : ds%dipoles(idir,:,ib_c,ib_v,ik,spin) = vv(:, idir)
456 : ! Hermitian conjugate
457 14700 : 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 30 : call vkbr_free(vkbr)
497 :
498 63 : if (nkpt < 1000 .or. (nkpt > 1000 .and. mod(ik, 200) == 0) .or. ik <= nproc) then
499 30 : write(msg,'(2(a,i0),a)')" k-point [", ik, "/", nkpt, "]"
500 30 : call cwtime_report(msg, cpu, wall, gflops)
501 : end if
502 :
503 : end do ! k-points
504 : end do ! spin
505 :
506 3 : call cwtime_report(msg, cpu_all, wall_all, gflops_all)
507 :
508 3 : ABI_FREE(ug_c)
509 3 : ABI_FREE(ug_v)
510 3 : ABI_FREE(kg_k)
511 3 : ABI_FREE(ihrc)
512 3 : ABI_FREE(cwaveprj)
513 3 : ABI_SFREE(distrib_mat)
514 3 : ABI_SFREE(distrib_diago)
515 :
516 3 : if (dtset%useria /= 666) then
517 3 : ABI_FREE(cg_c)
518 3 : ABI_FREE(cg_v)
519 3 : call ddkop%free()
520 : end if
521 :
522 : ! Gather the k-points computed by all processes
523 3 : call xmpi_sum_master(ds%dipoles, master, comm, ierr)
524 :
525 3 : if (ds%only_diago) then
526 0 : call xmpi_sum_master(ds%vdiago, master, comm, ierr)
527 : else
528 3 : 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 3 : 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 0 : call hdr%copy(tmp_hdr)
537 0 : 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 0 : do ii=1,3
563 0 : fname = strcat(prefix, '_', itoa(ii), "_EVK.nc")
564 0 : call wrtout(ab_out, sjoin("- Writing EVK file: ", fname, "for reduced direction:", itoa(ii)))
565 0 : NCF_CHECK_MSG(nctk_open_create(ncid, fname, xmpi_comm_self), "Creating EVK.nc file")
566 0 : tmp_hdr%pertcase = 3 * cryst%natom + ii
567 0 : NCF_CHECK(tmp_hdr%ncwrite(ncid, 43, nc_define=.True.))
568 0 : NCF_CHECK(cryst%ncwrite(ncid))
569 0 : NCF_CHECK(ebands%ncwrite(ncid))
570 : ncerr = nctk_def_arrays(ncid, [ &
571 : nctkarr_t('h1_matrix_elements', "dp", &
572 0 : "two, max_number_of_states, max_number_of_states, number_of_kpoints, number_of_spins")], defmode=.True.)
573 0 : NCF_CHECK(ncerr)
574 0 : NCF_CHECK(nctk_set_datamode(ncid))
575 0 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "h1_matrix_elements"), ds%dipoles(ii,:,:,:,:,:)))
576 0 : NCF_CHECK(nf90_close(ncid))
577 : end do
578 0 : call tmp_hdr%free()
579 : end if
580 :
581 3 : if (my_rank == master .and. dtset%prtvol > 0) then
582 3 : write(ab_out, "(2a)")ch10,"Writing velocity matrix elements (only diagonal terms, real part) for testing purpose:"
583 6 : do spin=1,nsppol
584 18 : do ik=1,min(nkpt, 4)
585 12 : write(ab_out, "(2(a,1x,i0),2x,2a)")"For spin: ", spin, ", ikbz: ", ik, ", kpt: ", trim(ktoa(wfd%kibz(:,ik)))
586 120 : do ib_c=bmin,min(bmin+8, bmax)
587 120 : write(ab_out, "(3(es16.6,2x))") ds%dipoles(:,1,ib_c,ib_c,ik,spin)
588 : end do
589 15 : 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 3 : call wfd%free()
599 3 : call ebands%free()
600 3 : call cryst%free()
601 3 : call hdr%free()
602 :
603 : ! Block all procs here so that we know output files are available when code returns.
604 3 : call xmpi_barrier(comm)
605 :
606 39 : 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 3 : subroutine ddkstore_free(self)
623 :
624 : !Arguments ------------------------------------
625 : class(ddkstore_t),intent(inout) :: self
626 : !************************************************************************
627 :
628 3 : ABI_SFREE(self%vdiago)
629 3 : ABI_SFREE(self%vmat)
630 3 : ABI_SFREE(self%dipoles)
631 :
632 3 : 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 74 : 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 74 : ABI_CHECK(dtset%usepaw == 0, "PAW not tested/implemented!")
712 :
713 74 : ddkop%inclvkb = dtset%inclvkb
714 74 : ddkop%usepaw = dtset%usepaw
715 74 : ddkop%ipert = cryst%natom + 1
716 74 : ddkop%dfpt_sciss = dtset%dfpt_sciss
717 74 : ddkop%mpw = mpw
718 :
719 962 : ddkop%rprimd = cryst%rprimd
720 74 : ddkop%mpi_enreg => mpi_enreg
721 :
722 : ! Not used because vlocal1 is not applied.
723 296 : nfft = product(ngfft(1:3))
724 296 : mgfft = maxval(ngfft(1:3))
725 :
726 296 : 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 222 : )
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 296 : call ddkop%rf_hamkq(idir)%init(cplex1, ddkop%gs_hamkq(idir), ddkop%ipert, has_e1kbsc=.true.)
742 : end do
743 :
744 74 : 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 361 : 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 361 : type(mpi_type) :: mpienreg_seq
785 : !arrays
786 : integer :: npwarr(nkpt1), dummy_nband(nkpt1*nsppol1)
787 : integer :: idir, nkpg, nkpg1, useylmgr1, optder !, nylmgr1
788 361 : real(dp),allocatable :: ylm_k(:,:),ylmgr1_k(:,:,:)
789 : !************************************************************************
790 :
791 361 : ABI_CHECK_ILEQ(npw_k, ddkop%mpw, "npw_k > mpw!")
792 1444 : ddkop%kpoint = kpoint
793 :
794 : ! Set up the spherical harmonics (Ylm) at k+q if useylm = 1
795 361 : useylmgr1 = 0; optder = 0
796 361 : if (psps%useylm == 1) then
797 18 : useylmgr1 = 1; optder = 1
798 : end if
799 :
800 2027986 : ABI_RECALLOC(ddkop%gh1c, (2, npw_k*dtset%nspinor, 3))
801 2027264 : ABI_RECALLOC(ddkop%gs1c, (2, npw_k*dtset%nspinor, 3))
802 :
803 1444 : ABI_MALLOC(ylm_k, (npw_k, psps%mpsang**2 * psps%useylm))
804 1805 : ABI_MALLOC(ylmgr1_k, (npw_k, 3+6*(optder/2), psps%mpsang**2*psps%useylm*useylmgr1))
805 :
806 361 : 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 1444 : do idir=1,3
816 1083 : call ddkop%htg(idir)%free()
817 :
818 : ! Continue to initialize the Hamiltonian
819 1083 : call ddkop%gs_hamkq(idir)%load_spin(spin, with_nonlocal=.true.)
820 1083 : 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 1444 : ddkop%htg(idir)%ph3d, ddkop%htg(idir)%ph3d1) ! Out
831 : end do
832 :
833 361 : ABI_FREE(ylm_k)
834 361 : ABI_FREE(ylmgr1_k)
835 :
836 361 : 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 1688 : 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 1688 : ddkop%eig0nk = eig0nk
884 :
885 1688 : 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 1588 : optnl = 2; if (ddkop%inclvkb == 0) optnl = 0
891 :
892 1588 : eshift = ddkop%eig0nk - ddkop%dfpt_sciss
893 6352 : do idir=1,3
894 4764 : 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 11116 : 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 1688 : 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 6076 : 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 3038 : if (ddkop%usepaw == 0) then
962 : ! <u_(iband,k+q)^(0)|H_(k+q,k)^(1)|u_(jband,k)^(0)> (NC psps)
963 12152 : do idir=1,3
964 9114 : dotarr = cg_zdotc(npw_k * nspinor, brag, ddkop%gh1c(:,:,idir))
965 9114 : if (istwf_k > 1) then
966 : !dum = two * j_dpc * AIMAG(dum); if (vkbr%istwfk==2) dum = dum - j_dpc * AIMAG(gamma_term)
967 0 : doti = two * dotarr(2)
968 0 : 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 0 : doti = doti - (brag(1) * ddkop%gh1c(2,1,idir) - brag(2) * ddkop%gh1c(1,1,idir))
972 : end if
973 0 : dotarr(2) = doti; dotarr(1) = zero
974 : end if
975 30380 : 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 3038 : 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 3038 : ABI_ERROR(sjoin("Invalid vaue for mode:", mode))
993 : end select
994 :
995 3038 : 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 74 : 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 74 : ABI_SFREE(ddkop%gh1c)
1177 74 : ABI_SFREE(ddkop%gs1c)
1178 :
1179 296 : do idir=1,3
1180 222 : call ddkop%gs_hamkq(idir)%free()
1181 222 : call ddkop%htg(idir)%free()
1182 296 : call ddkop%rf_hamkq(idir)%free()
1183 : end do
1184 :
1185 74 : ddkop%mpi_enreg => null()
1186 :
1187 74 : 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 1305 : subroutine ham_targets_free(ddkop)
1202 :
1203 : !Arguments ------------------------------------
1204 : class(ham_targets_t),intent(inout) :: ddkop
1205 : !************************************************************************
1206 :
1207 1305 : ABI_SFREE(ddkop%ffnlk)
1208 1305 : ABI_SFREE(ddkop%ffnl1)
1209 1305 : ABI_SFREE(ddkop%kpg_k)
1210 1305 : ABI_SFREE(ddkop%kpg1_k)
1211 1305 : ABI_SFREE(ddkop%dkinpw)
1212 1305 : ABI_SFREE(ddkop%kinpw1)
1213 1305 : ABI_SFREE(ddkop%ph3d)
1214 1305 : ABI_SFREE(ddkop%ph3d1)
1215 :
1216 1305 : end subroutine ham_targets_free
1217 : !!***
1218 :
1219 : !----------------------------------------------------------------------
1220 :
1221 5476 : end module m_ddk
1222 : !!***
|
