Line data Source code
1 : !!****m* ABINIT/m_skw
2 : !! NAME
3 : !! m_skw
4 : !!
5 : !! FUNCTION
6 : !! Shankland-Koelling-Wood Fourier interpolation scheme.
7 : !!
8 : !! COPYRIGHT
9 : !! Copyright (C) 2008-2026 ABINIT group (MG)
10 : !! This file is distributed under the terms of the
11 : !! GNU General Public License, see ~abinit/COPYING
12 : !! or http://www.gnu.org/copyleft/gpl.txt .
13 : !!
14 : !! SOURCE
15 :
16 : #if defined HAVE_CONFIG_H
17 : #include "config.h"
18 : #endif
19 :
20 : #include "abi_common.h"
21 :
22 : module m_skw
23 :
24 : use defs_basis
25 : use m_errors
26 : use m_abicore
27 : use m_xmpi
28 : use m_crystal
29 : use m_sort
30 : use m_nctk
31 : use netcdf
32 :
33 : use m_fstrings, only : itoa, sjoin, ktoa, yesno, ftoa
34 : use m_special_funcs, only : abi_derfc
35 : use m_time, only : cwtime, cwtime_report
36 : use m_numeric_tools, only : imax_loc
37 : use m_bz_mesh, only : isamek
38 : use m_gsphere, only : get_irredg
39 :
40 : implicit none
41 :
42 : private
43 : !!***
44 :
45 : !----------------------------------------------------------------------
46 :
47 : !!****t* m_skw/skw_t
48 : !! NAME
49 : !! skw_t
50 : !!
51 : !! FUNCTION
52 : !! Object implementing the Shankland-Koelling-Wood Fourier interpolation scheme.
53 : !! It can be used to interpolate functions in k-space with the periodicity of the
54 : !! reciprocal lattice and satisfying F(k) = F(Sk) for each rotation S
55 : !! belonging to the point group of the crystal. For readability reason,
56 : !! the names of the variables are chosen assuming we are interpolating electronic eigenvalues
57 : !! but the same object can be used to interpolate phonons as well. Just use nsppol=1 and nband = 3 * natom
58 : !!
59 : !! SOURCE
60 :
61 : type,public :: skw_t
62 :
63 : integer :: cplex
64 : ! 1 if time-reversal symmetry can be used, 2 otherwise.
65 :
66 : integer :: nr
67 : ! Number of star functions.
68 :
69 : integer :: nkpt
70 : ! Number of ab-initio k-points.
71 :
72 : integer :: ptg_nsym
73 : ! Number of operations in the point group.
74 :
75 : logical :: has_inversion
76 : ! True if the point group contains spatial inversion.
77 :
78 : integer :: band_block(2)
79 : ! Initial and final band index.
80 :
81 : integer :: bcount
82 : ! Number of bands
83 :
84 : integer :: nsppol
85 : ! Number of independent spin polarizations.
86 :
87 : integer,allocatable :: rpts(:,:)
88 : ! rpts(3, nr)
89 : ! Real-space lattice points (in reduced coordinates) ordered with non-decreasing length.
90 :
91 : integer,allocatable :: ptg_symrel(:,:,:)
92 : ! ptg_symrel(3,3,ptg_nsym)
93 : ! operations of the point group (real space).
94 :
95 : integer,allocatable :: ptg_symrec(:,:,:)
96 : ! ptg_symrec(3,3,ptg_nsym)
97 : ! operations of the point group (reciprocal space).
98 :
99 : complex(dp),allocatable :: coefs(:,:,:)
100 : ! coefs(nr, bcount, nsppol).
101 :
102 : complex(dp),allocatable :: cached_srk(:)
103 : ! cached_srk(%nr)
104 : ! The star function for cached_kpt (used in skw_eval_bks).
105 : real(dp) :: cached_kpt(3)
106 :
107 : complex(dp),allocatable :: cached_srk_dk1(:,:)
108 : ! cached_srk_dk1(%nr, 3)
109 : ! The 1d derivative wrt k of the star function for cached_kpt_dk1 (used in skw_eval_bks).
110 : real(dp) :: cached_kpt_dk1(3)
111 :
112 : complex(dp),allocatable :: cached_srk_dk2(:,:,:)
113 : ! cached_srk_dk2(%nr,3,3)
114 : ! The 2d derivatives wrt k of the star function for cached_kpt_dk2 (used in skw_eval_bks).
115 : real(dp) :: cached_kpt_dk2(3)
116 :
117 : contains
118 :
119 : procedure :: print => skw_print
120 : ! Print info about object.
121 :
122 : procedure :: ncwrite => skw_ncwrite
123 : ! Write the object in netcdf format
124 :
125 : procedure :: eval_bks => skw_eval_bks
126 : ! Interpolate eigenvalues, 1st, 2nd derivates wrt k, at an arbitrary k-point.
127 :
128 : procedure :: free => skw_free
129 : ! Free memory.
130 :
131 : end type skw_t
132 : !!***
133 :
134 : public :: skw_new ! Create new object.
135 :
136 : CONTAINS !=====================================================================================
137 : !!***
138 :
139 : !!****f* m_skw/skw_new
140 : !! NAME
141 : !! skw_new
142 : !!
143 : !! FUNCTION
144 : !! Initialize the object.
145 : !!
146 : !! INPUTS
147 : !! cryst<crystal_t>=Crystalline structure.
148 : !! params(:)
149 : !! params(1): Ratio between star functions and ab-initio k-points.
150 : !! params(2:3): Activate Fourier filtering (Eq 9 of PhysRevB.61.1639 [[cite:Uehara2000]]) if params(2) > tol6
151 : !! params(2)=rcut, params(3) = rsigma
152 : !! cplex=1 if time reversal can be used, 2 otherwise.
153 : !! nband=Total Number of bands in the eig array.
154 : !! nkpt=Number of ab-initio k-points.
155 : !! nsppol=Number of independent spin polarizations.
156 : !! kpts(3,nkpt)=ab-initio k-points in reduced coordinates.
157 : !! eig(nband,nkpt,nsppol)=ab-initio eigenvalues.
158 : !! band_block(2)=Initial and final band index to interpolate. If [0,0], all bands are used
159 : !! This is a global variable i.e. all MPI procs MUST call the routine with the same value.
160 : !! comm=MPI communicator
161 : !!
162 : !! SOURCE
163 :
164 12 : type(skw_t) function skw_new(cryst, params, cplex, nband, nkpt, nsppol, kpts, eig, band_block, comm) result(new)
165 :
166 : !Arguments ------------------------------------
167 : !scalars
168 : integer,intent(in) :: cplex,nband,nkpt,nsppol,comm
169 : real(dp),intent(in) :: params(:)
170 : type(crystal_t),intent(in) :: cryst
171 : !arrays
172 : integer,intent(in) :: band_block(2)
173 : real(dp),intent(in) :: kpts(3,nkpt)
174 : real(dp),intent(in) :: eig(nband,nkpt,nsppol)
175 :
176 : !Local variables-------------------------------
177 : !scalars
178 : integer,parameter :: master=0,prtvol=1
179 : integer :: my_rank,nprocs,cnt,bstop,bstart,bcount,lwork
180 : integer :: ir,ik,ib,ii,jj,nr,band,spin,ierr,lpratio,nrwant
181 : real(dp),parameter :: c1=0.25_dp,c2=0.25_dp
182 : real(dp) :: r2,r2min,mare,mae_meV,adiff_meV,rel_err,rcut,rsigma
183 : real(dp) :: cpu_tot,wall_tot,gflops_tot,cpu,wall,gflops,rval
184 : character(len=500) :: fmt,msg
185 : !arrays
186 : integer :: rmax(3)
187 6 : integer,allocatable :: ipiv(:)
188 : real(dp) :: list2(2)
189 6 : real(dp),allocatable :: r2vals(:),inv_rhor(:),oeig(:)
190 6 : complex(dp),allocatable :: srk(:,:),hmat(:,:),lambda(:,:,:),work(:)
191 : ! *********************************************************************
192 :
193 6 : ABI_CHECK(nkpt > 1, sjoin("nkpt must be > 1 but got:", itoa(nkpt)))
194 :
195 6 : call cwtime(cpu_tot, wall_tot, gflops_tot, "start")
196 :
197 6 : nprocs = xmpi_comm_size(comm); my_rank = xmpi_comm_rank(comm)
198 :
199 : ! Get slice of bands to be treated.
200 18 : new%band_block = band_block; if (all(band_block == 0)) new%band_block = [1, nband]
201 6 : bstart = new%band_block(1); bstop = new%band_block(2); bcount = bstop - bstart + 1
202 6 : new%cplex = cplex; new%nkpt = nkpt; new%nsppol = nsppol; new%bcount = bcount
203 :
204 : ! Get point group operations.
205 6 : call cryst%get_point_group(new%ptg_nsym, new%ptg_symrel, new%ptg_symrec, new%has_inversion, include_timrev=cplex==1)
206 :
207 : ! -----------------------
208 : ! Find nrwant star points
209 : ! -----------------------
210 6 : lpratio = int(abs(params(1)))
211 6 : ABI_CHECK(lpratio > 0, "lpratio must be > 0")
212 24 : rmax = nint((one + (lpratio * new%nkpt * new%ptg_nsym) / two) ** third)
213 6 : if (new%has_inversion) then
214 24 : rmax = nint((one + (lpratio * new%nkpt * new%ptg_nsym / 2) / two) ** third)
215 : end if
216 6 : nrwant = lpratio * new%nkpt
217 :
218 6 : call cwtime(cpu, wall, gflops, "start")
219 0 : do
220 6 : call find_rstar_gen(new, cryst, nrwant, rmax, r2vals, comm)
221 6 : if (new%nr >= nrwant) then
222 : !write(std_out,*)"Entered with rmax", rmax," abs(skw%rpts(last)): ", abs(new%rpts(:,new%nr))
223 : exit
224 : end if
225 0 : write(std_out,*)"rmax: ", rmax," was not large enough to find ", nrwant," R-star points."
226 0 : rmax = 2 * rmax
227 0 : write(std_out,*)"Will try again with enlarged rmax: ",rmax
228 0 : ABI_FREE(r2vals)
229 : end do
230 6 : nr = new%nr
231 6 : call cwtime_report(" find_rstar_gen", cpu, wall, gflops)
232 :
233 6 : if (my_rank == master) call new%print(std_out)
234 :
235 : ! Compute (inverse) roughness function.
236 6 : r2min = r2vals(2)
237 18 : ABI_MALLOC(inv_rhor, (nr))
238 341 : do ir=1,nr
239 335 : r2 = r2vals(ir)
240 341 : inv_rhor(ir) = one / ((one - c1 * r2/r2min)**2 + c2 * (r2 / r2min)**3)
241 : ! TODO: Test the two versions.
242 : !if (params(1) < zero) inv_rhor(ir) = one / (c1 * r2 + c2 * r2**2)
243 : end do
244 :
245 : ! Construct star functions for the ab-initio k-points.
246 24 : ABI_MALLOC(srk, (nr, nkpt))
247 73 : do ik=1,nkpt
248 73 : call mkstar(new, kpts(:,ik), srk(:,ik))
249 : end do
250 :
251 : ! Build H(k,k') matrix (Hermitian)
252 1090 : ABI_CALLOC(hmat, (nkpt-1, nkpt-1))
253 : cnt = 0
254 67 : do jj=1,nkpt-1
255 600 : do ii=1,jj
256 : !do ii=1,nkpt-1
257 533 : cnt = cnt + 1; if (mod(cnt, nprocs) /= my_rank) cycle ! mpi parallelism.
258 63861 : do ir=2,nr
259 : hmat(ii, jj) = hmat(ii, jj) + &
260 63800 : (srk(ir, ii) - srk(ir, nkpt)) * conjg(srk(ir, jj) - srk(ir, nkpt)) * inv_rhor(ir)
261 : end do
262 : end do
263 : end do
264 6 : call xmpi_sum(hmat, comm, ierr)
265 :
266 30 : ABI_MALLOC(lambda, (nkpt-1, bcount, nsppol))
267 12 : do spin=1,nsppol
268 88 : do ib=1,bcount
269 76 : band = ib + bstart - 1
270 754 : lambda(:,ib,spin) = eig(band,1:nkpt-1,spin) - eig(band,nkpt,spin)
271 : end do
272 : end do
273 :
274 : ! Solve all bands and spins at once [[cite:Pickett1988]]
275 6 : call wrtout(std_out, " Solving system of linear equations to get lambda coeffients (eq. 10 of PRB 38 2721)...", do_flush=.True.)
276 6 : call cwtime(cpu, wall, gflops, "start")
277 18 : ABI_MALLOC(ipiv, (nkpt-1))
278 :
279 : if (.False.) then
280 : ! General complex.
281 : call zgesv(nkpt-1, bcount*nsppol, hmat, nkpt-1, ipiv, lambda, nkpt-1, ierr)
282 : ABI_CHECK(ierr == 0, sjoin("ZGESV returned:", itoa(ierr)))
283 : else
284 : ! Hermitian version
285 67 : do ii=1,nkpt-1
286 67 : hmat(ii, ii) = real(hmat(ii, ii))
287 : end do
288 6 : lwork = -1
289 6 : ABI_MALLOC(work, (1))
290 6 : call zhesv("U", nkpt-1, bcount*nsppol, hmat, nkpt-1, ipiv, lambda, nkpt-1, work, lwork, ierr)
291 6 : lwork = nint(real(work(1)))
292 6 : ABI_FREE(work)
293 18 : ABI_MALLOC(work, (lwork))
294 6 : call zhesv("U", nkpt-1, bcount*nsppol, hmat, nkpt-1, ipiv, lambda, nkpt-1, work, lwork, ierr)
295 6 : ABI_CHECK(ierr == 0, sjoin("ZHESV returned:", itoa(ierr)))
296 6 : ABI_FREE(work)
297 : end if
298 6 : call cwtime_report(" ZHESV", cpu, wall, gflops)
299 :
300 : ! Compute coefficients
301 30 : ABI_MALLOC(new%coefs, (nr,bcount,nsppol))
302 :
303 12 : do spin=1,nsppol
304 88 : do ib=1,bcount
305 76 : band = ib + bstart - 1
306 3740 : do ir=2,nr
307 52768 : new%coefs(ir,ib,spin) = inv_rhor(ir) * dot_product(srk(ir,:nkpt-1) - srk(ir,nkpt), lambda(:nkpt-1, ib, spin))
308 : !new%coefs(ir,ib,spin) = inv_rhor(ir) * dot_product(lambda(:nkpt-1, ib, spin), conjg(srk(ir,:) - srk(ir,nkpt)))
309 : !new%coefs(ir,ib,spin) = inv_rhor(ir) * dot_product(lambda(:nkpt-1, ib, spin), conjg(srk(ir,:) - srk(ir,1)))
310 : end do
311 3746 : new%coefs(1,ib,spin) = eig(band,nkpt,spin) - dot_product(conjg(new%coefs(2:nr, ib,spin)), srk(2:nr, nkpt))
312 : end do
313 : end do
314 :
315 : ! Filter high-frequency.
316 6 : if (params(2) > tol6) then
317 0 : rcut = params(2) * sqrt(r2vals(new%nr))
318 0 : rsigma = params(3); if (rsigma <= zero) rsigma = five
319 0 : call wrtout(std_out," Applying filter (Eq 9 of PhysRevB.61.1639)") ! [[cite:Uehara2000]]
320 0 : do ir=2,nr
321 0 : new%coefs(ir,:,:) = new%coefs(ir,:,:) * half * abi_derfc((sqrt(r2vals(ir)) - rcut) / rsigma)
322 : end do
323 : end if
324 :
325 : ! Prepare workspace arrays for star functions.
326 24 : new%cached_kpt = huge(one)
327 18 : ABI_MALLOC(new%cached_srk, (new%nr))
328 24 : new%cached_kpt_dk1 = huge(one)
329 18 : ABI_MALLOC(new%cached_srk_dk1, (new%nr, 3))
330 24 : new%cached_kpt_dk2 = huge(one)
331 24 : ABI_MALLOC(new%cached_srk_dk2, (new%nr, 3, 3))
332 :
333 6 : ABI_FREE(r2vals)
334 6 : ABI_FREE(srk)
335 6 : ABI_FREE(inv_rhor)
336 6 : ABI_FREE(hmat)
337 6 : ABI_FREE(lambda)
338 6 : ABI_FREE(ipiv)
339 :
340 : ! Compare ab-initio data with interpolated results.
341 18 : ABI_MALLOC(oeig, (bcount))
342 6 : fmt = sjoin("(a,", itoa(bcount), "(es12.4))")
343 6 : bstop = bstart + bcount - 1
344 6 : mare = zero; mae_meV = zero; cnt = 0
345 6 : call wrtout(std_out, ch10//" Comparing ab-initio energies with SKW interpolated results...")
346 12 : do spin=1,nsppol
347 79 : do ik=1,nkpt
348 67 : cnt = cnt + 1; if (mod(cnt, nprocs) /= my_rank) cycle ! mpi parallelism.
349 :
350 815 : do ib=1,bcount
351 748 : band = ib + new%band_block(1) - 1
352 748 : call new%eval_bks(band, kpts(:,ik), spin, oeig(ib))
353 :
354 748 : adiff_meV = abs(eig(band,ik,spin) - oeig(ib)); rel_err = zero
355 748 : if (abs(eig(band,ik,spin)) > tol16) rel_err = adiff_meV / abs(eig(band,ik,spin))
356 748 : rel_err = 100 * rel_err; adiff_meV = adiff_meV * Ha_meV
357 815 : mae_meV = mae_meV + adiff_meV; mare = mare + rel_err
358 : end do
359 :
360 6 : if (prtvol > 0) then
361 815 : ib = imax_loc(eig(bstart:bstop,ik,spin) - oeig)
362 67 : rval = (eig(bstart+ib-1,ik,spin) - oeig(ib)) * Ha_meV
363 : write(std_out,"(a,es12.4,2a)") &
364 67 : " SKW maxerr: ", rval, &
365 134 : " (meV), kpt: ", sjoin(ktoa(kpts(:,ik)), "band:",itoa(bstart+ib-1),", spin: ", itoa(spin))
366 : !write(std_out,fmt)"-- ref ", eig(bstart:bstop,ik,spin) * Ha_meV
367 : !write(std_out,fmt)"-- int ", oeig * Ha_meV
368 : !call vdiff%eval(1, bcount, eig(bstart:bstop,ik,spin), oeig, one, unit=std_out))
369 : end if
370 : end do
371 : end do
372 6 : ABI_FREE(oeig)
373 :
374 : ! Issue warning if error too large.
375 18 : list2 = [mare, mae_meV]; call xmpi_sum(list2, comm, ierr); mare = list2(1); mae_meV = list2(2)
376 6 : cnt = bcount * nkpt * nsppol; mare = mare / cnt; mae_meV = mae_meV / cnt
377 6 : write(std_out,"(2(a,es12.4),a,/)")" MARE: ",mare, ", MAE: ", mae_meV, " (meV)"
378 6 : if (mae_meV > ten) then
379 : write(msg,"(2a,2(a,es12.4),a)") &
380 0 : "Large error in SKW interpolation!",ch10," MARE: ",mare, ", MAE: ", mae_meV, " (meV)"
381 0 : call wrtout(ab_out, msg)
382 0 : ABI_WARNING(msg)
383 : end if
384 :
385 6 : call cwtime_report(" skw_new", cpu_tot, wall_tot, gflops_tot, end_str=ch10)
386 :
387 48 : end function skw_new
388 : !!***
389 :
390 : !----------------------------------------------------------------------
391 :
392 : !!****f* m_skw/skw_print
393 : !! NAME
394 : !! skw_print
395 : !!
396 : !! FUNCTION
397 : !! Print info on object
398 : !!
399 : !! INPUTS
400 : !! unt=Fortran unit number.
401 : !!
402 : !! OUTPUT
403 : !! only writing
404 : !!
405 : !! SOURCE
406 :
407 6 : subroutine skw_print(skw, unt)
408 :
409 : !Arguments ------------------------------------
410 : !scalars
411 : class(skw_t),intent(in) :: skw
412 : integer,intent(in) :: unt
413 :
414 : ! *********************************************************************
415 :
416 6 : write(unt,"(a)")" === Shankland-Koelling-Wood Fourier interpolation scheme ==="
417 6 : write(unt,"(a)")sjoin(" nsppol", itoa(skw%nsppol), ", cplex:", itoa(skw%cplex))
418 6 : write(unt,"(a)")sjoin(" Number of ab-initio k-points:", itoa(skw%nkpt))
419 6 : write(unt,"(a)")sjoin(" Number of star functions:", itoa(skw%nr))
420 6 : write(unt,"(a)")sjoin(" Stars/Nk ratio:", ftoa(skw%nr * one / skw%nkpt))
421 6 : write(unt,"(a)")sjoin(" Has spatial inversion:", yesno(skw%has_inversion))
422 :
423 6 : end subroutine skw_print
424 : !!***
425 :
426 : !----------------------------------------------------------------------
427 :
428 : !!****f* m_skw/skw_ncwrite
429 : !! NAME
430 : !! skw_ncwrite
431 : !!
432 : !! FUNCTION
433 : !! Write the object in netcdf format
434 : !!
435 : !! INPUTS
436 : !! ncid=NC file handle.
437 : !! [prefix]=String prepended to netcdf dimensions/variables (HDF5 poor-man groups)
438 : !! "skw" if not specified.
439 : !!
440 : !! OUTPUT
441 : !! Only writing
442 : !!
443 : !! SOURCE
444 :
445 0 : integer function skw_ncwrite(self, ncid, prefix) result(ncerr)
446 :
447 : !Arguments ------------------------------------
448 : !scalars
449 : class(skw_t),intent(in) :: self
450 : integer,intent(in) :: ncid
451 : character(len=*),optional,intent(in) :: prefix
452 :
453 : !Local variables-------------------------------
454 : !scalars
455 : character(len=500) :: prefix_
456 : !arrays
457 0 : real(dp),allocatable :: real_coefs(:,:,:,:)
458 : ! *************************************************************************
459 :
460 0 : prefix_ = "skw"; if (present(prefix)) prefix_ = trim(prefix)
461 :
462 : ! Define dimensions.
463 : ncerr = nctk_def_dims(ncid, [ &
464 : nctkdim_t("nr", self%nr), nctkdim_t("nkpt", self%nkpt), nctkdim_t("bcount", self%bcount), &
465 : nctkdim_t("nsppol", self%nsppol)], &
466 0 : defmode=.True., prefix=prefix_)
467 0 : NCF_CHECK(ncerr)
468 :
469 : ncerr = nctk_def_arrays(ncid, [ &
470 : ! Atomic structure and symmetry operations
471 : nctkarr_t("rpts", "dp", "three, number_of_cartesian_directions, number_of_vectors"), &
472 : nctkarr_t("kpts", "dp", "three, nkpt"), &
473 : nctkarr_t("coefs", "dp", "two, nr, bcount, nsppol") &
474 0 : ], prefix=prefix_)
475 0 : NCF_CHECK(ncerr)
476 :
477 : ! Write data.
478 0 : NCF_CHECK(nctk_set_datamode(ncid))
479 0 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, pre("rpts")), self%rpts))
480 : !NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, pre("kpts")), self%kpts))
481 0 : ABI_MALLOC(real_coefs, (2, self%nr, self%bcount, self%nsppol))
482 0 : real_coefs(1,:,:,:) = real(self%coefs); real_coefs(2,:,:,:) = aimag(self%coefs)
483 0 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, pre("coefs")), real_coefs))
484 0 : ABI_FREE(real_coefs)
485 :
486 : contains
487 0 : pure function pre(istr) result(ostr)
488 : character(len=*),intent(in) :: istr
489 : character(len=len_trim(prefix_) + len_trim(istr)+1) :: ostr
490 0 : ostr = trim(prefix_) // trim(istr)
491 0 : end function pre
492 :
493 : end function skw_ncwrite
494 : !!***
495 :
496 : !!****f* m_skw/skw_eval_bks
497 : !! NAME
498 : !! skw_eval_bks
499 : !!
500 : !! FUNCTION
501 : !! Interpolate the energies for an arbitrary k-point and spin with slow FT.
502 : !!
503 : !! INPUTS
504 : !! band=Band index (global index associated to the input eigenvalues, i.e. independent of band_block)
505 : !! kpt(3)=K-point in reduced coordinates.
506 : !! spin=Spin index.
507 : !!
508 : !! OUTPUT
509 : !! oeig=interpolated eigenvalues
510 : !! Note that oeig is not necessarily sorted in ascending order.
511 : !! The routine does not reorder the interpolated eigenvalues
512 : !! to be consistent with the interpolation of the derivatives.
513 : !! [oder1(3)]=First-order derivatives wrt k in reduced coordinates.
514 : !! [oder2(3,3)]=Second-order derivatives wrt k in reduced coordinates.
515 : !!
516 : !! SOURCE
517 :
518 13434 : subroutine skw_eval_bks(skw, band, kpt, spin, oeig, oder1, oder2)
519 :
520 : !Arguments ------------------------------------
521 : !scalars
522 : integer,intent(in) :: band,spin
523 : class(skw_t),intent(inout) :: skw
524 : !arrays
525 : real(dp),intent(in) :: kpt(3)
526 : real(dp),intent(out) :: oeig
527 : real(dp),optional,intent(out) :: oder1(3),oder2(3,3)
528 :
529 : !Local variables-------------------------------
530 : !scalars
531 : integer :: ii,jj,ib
532 : ! *********************************************************************
533 :
534 13434 : ib = band - skw%band_block(1) + 1
535 13434 : ABI_CHECK(ib >= 1 .and. ib <= skw%bcount, sjoin("out of range band:", itoa(band)))
536 :
537 : ! Compute star function for this k-point (if not already in memory)
538 50579 : if (any(kpt /= skw%cached_kpt)) then
539 1098 : call mkstar(skw, kpt, skw%cached_srk)
540 4392 : skw%cached_kpt = kpt
541 : end if
542 :
543 854314 : oeig = dot_product(conjg(skw%coefs(:,ib,spin)), skw%cached_srk)
544 :
545 : ! TODO: Test Derivatives
546 13434 : if (present(oder1)) then
547 : ! Compute first-order derivatives.
548 0 : if (any(kpt /= skw%cached_kpt_dk1)) then
549 0 : call mkstar_dk1(skw, kpt, skw%cached_srk_dk1)
550 0 : skw%cached_kpt_dk1 = kpt
551 : end if
552 :
553 0 : do ii=1,3
554 0 : oder1(ii) = dot_product(conjg(skw%coefs(:,ib,spin)), skw%cached_srk_dk1(:,ii)) * two_pi
555 : end do
556 : end if
557 :
558 13434 : if (present(oder2)) then
559 : ! Compute second-order derivatives.
560 0 : if (any(kpt /= skw%cached_kpt_dk2)) then
561 0 : call mkstar_dk2(skw, kpt, skw%cached_srk_dk2)
562 0 : skw%cached_kpt_dk2 = kpt
563 : end if
564 :
565 0 : oder2 = zero
566 0 : do jj=1,3
567 0 : do ii=1,jj
568 0 : oder2(ii, jj) = dot_product(conjg(skw%coefs(:,ib,spin)), skw%cached_srk_dk2(:,ii,jj)) * two_pi**2
569 0 : if (ii /= jj) oder2(jj, ii) = oder2(ii, jj)
570 : end do
571 : end do
572 : end if
573 :
574 13434 : end subroutine skw_eval_bks
575 : !!***
576 :
577 : !----------------------------------------------------------------------
578 :
579 : !!****f* m_skw/skw_free
580 : !! NAME
581 : !! skw_free
582 : !!
583 : !! FUNCTION
584 : !! Free memory
585 : !!
586 : !! SOURCE
587 :
588 6 : subroutine skw_free(skw)
589 :
590 : !Arguments ------------------------------------
591 : !scalars
592 : class(skw_t),intent(inout) :: skw
593 :
594 : ! *********************************************************************
595 :
596 6 : ABI_SFREE(skw%rpts)
597 6 : ABI_SFREE(skw%ptg_symrel)
598 6 : ABI_SFREE(skw%ptg_symrec)
599 6 : ABI_SFREE(skw%coefs)
600 :
601 6 : ABI_SFREE(skw%cached_srk)
602 24 : skw%cached_kpt = huge(one)
603 6 : ABI_SFREE(skw%cached_srk_dk1)
604 24 : skw%cached_kpt_dk1 = huge(one)
605 6 : ABI_SFREE(skw%cached_srk_dk2)
606 24 : skw%cached_kpt_dk2 = huge(one)
607 :
608 6 : end subroutine skw_free
609 : !!***
610 :
611 : !----------------------------------------------------------------------
612 :
613 : !!****f* m_skw/mkstar
614 : !! NAME
615 : !! mkstar
616 : !!
617 : !! FUNCTION
618 : !! Compute the star function for k-point kpt
619 : !!
620 : !! INPUTS
621 : !! kpt(3)=K-point in reduced coordinates.
622 : !!
623 : !! OUTPUT
624 : !! srk(%nr)=Star function for this k-point.
625 : !!
626 : !! SOURCE
627 :
628 1165 : subroutine mkstar(skw, kpt, srk)
629 :
630 : !Arguments ------------------------------------
631 : !scalars
632 : type(skw_t),intent(in) :: skw
633 : !arrays
634 : real(dp),intent(in) :: kpt(3)
635 : complex(dp),intent(out) :: srk(skw%nr)
636 :
637 : !Local variables-------------------------------
638 : !scalars
639 : integer :: ir,isym
640 : !arrays
641 : real(dp) :: sk(3)
642 : ! *********************************************************************
643 :
644 91745 : srk = zero
645 57085 : do isym=1,skw%ptg_nsym
646 1621680 : sk = two_pi * matmul(transpose(skw%ptg_symrel(:,:,isym)), kpt)
647 4404925 : do ir=1,skw%nr
648 17447280 : srk(ir) = srk(ir) + exp(j_dpc * dot_product(sk, skw%rpts(:,ir)))
649 : end do
650 : end do
651 91745 : srk = srk / skw%ptg_nsym
652 :
653 1165 : end subroutine mkstar
654 : !!***
655 :
656 : !----------------------------------------------------------------------
657 :
658 : !!****f* m_skw/mkstar_dk1
659 : !! NAME
660 : !! mkstar_dk1
661 : !!
662 : !! FUNCTION
663 : !! Compute the 1st derivative of the star function wrt k
664 : !!
665 : !! INPUTS
666 : !! kpt(3)=K-point in reduced coordinates.
667 : !!
668 : !! OUTPUT
669 : !! srk_dk1(%nr,3)=Derivative of the star function wrt k in reduced coordinates.
670 : !!
671 : !! SOURCE
672 :
673 0 : subroutine mkstar_dk1(skw, kpt, srk_dk1)
674 :
675 : !Arguments ------------------------------------
676 : !scalars
677 : type(skw_t),intent(in) :: skw
678 : !arrays
679 : real(dp),intent(in) :: kpt(3)
680 : complex(dp),intent(out) :: srk_dk1(skw%nr,3)
681 :
682 : !Local variables-------------------------------
683 : !scalars
684 : integer :: ir,isym
685 : !arrays
686 : real(dp) :: sk(3)
687 0 : complex(dp) :: work(3,skw%nr)
688 : ! *********************************************************************
689 :
690 0 : work = zero
691 0 : do isym=1,skw%ptg_nsym
692 0 : sk = two_pi * matmul(transpose(skw%ptg_symrel(:,:,isym)), kpt)
693 0 : do ir=1,skw%nr
694 : work(:,ir) = work(:,ir) + exp(j_dpc * dot_product(sk, skw%rpts(:,ir))) * &
695 0 : matmul(skw%ptg_symrel(:,:,isym), skw%rpts(:,ir))
696 : end do
697 : end do
698 0 : work = j_dpc * work / skw%ptg_nsym
699 0 : srk_dk1 = transpose(work)
700 :
701 0 : end subroutine mkstar_dk1
702 : !!***
703 :
704 : !----------------------------------------------------------------------
705 :
706 : !!****f* m_skw/mkstar_dk2
707 : !! NAME
708 : !! mkstar_dk2
709 : !!
710 : !! FUNCTION
711 : !! Compute the 2st derivatives of the star function wrt k
712 : !!
713 : !! INPUTS
714 : !! kpt(3)=K-point in reduced coordinates.
715 : !!
716 : !! OUTPUT
717 : !! srk_dk2(%nr,3,3)=2nd derivatives of the star function wrt k in reduced coordinates.
718 : !!
719 : !! SOURCE
720 :
721 0 : subroutine mkstar_dk2(skw, kpt, srk_dk2)
722 :
723 : !Arguments ------------------------------------
724 : !scalars
725 : type(skw_t),intent(in) :: skw
726 : !arrays
727 : real(dp),intent(in) :: kpt(3)
728 : complex(dp),intent(out) :: srk_dk2(skw%nr,3,3)
729 :
730 : !Local variables-------------------------------
731 : !scalars
732 : integer :: ir,isym,ii,jj
733 : complex(dp) :: eiskr
734 : !arrays
735 : integer :: sr(3)
736 : real(dp) :: sk(3)
737 0 : complex(dp) :: work(3,3,skw%nr)
738 : ! *********************************************************************
739 :
740 0 : work = zero
741 0 : do isym=1,skw%ptg_nsym
742 0 : sk = two_pi * matmul(transpose(skw%ptg_symrel(:,:,isym)), kpt)
743 0 : do ir=1,skw%nr
744 0 : sr = matmul(skw%ptg_symrel(:,:,isym), skw%rpts(:,ir))
745 0 : eiskr = exp(j_dpc * dot_product(sk, skw%rpts(:,ir)))
746 0 : do jj=1,3
747 0 : do ii=1,jj
748 0 : work(ii,jj,ir) = work(ii,jj,ir) + eiskr * sr(ii) * sr(jj)
749 : end do
750 : end do
751 : end do
752 : end do
753 0 : work = - work / skw%ptg_nsym
754 :
755 0 : do jj=1,3
756 0 : do ii=1,jj
757 0 : srk_dk2(:, ii, jj) = work(ii, jj, :)
758 0 : if (ii /= jj) srk_dk2(:,jj,ii) = work(:,ii,jj)
759 : end do
760 : end do
761 :
762 0 : end subroutine mkstar_dk2
763 : !!***
764 :
765 : !----------------------------------------------------------------------
766 :
767 : !!****f* m_skw/find_rstar_gen
768 : !! NAME
769 : !! find_rstar_gen
770 : !!
771 : !! FUNCTION
772 : !! Find the R-space points generating the stars.
773 : !! Set skw%nr and skw%rpts.
774 : !!
775 : !! INPUTS
776 : !! cryst<crystal_t>=Crystalline structure.
777 : !! nrwant=Number of R-space points wanted
778 : !! rmax(3)=Max reduced components of supercell.
779 : !! comm=MPI communicator.
780 : !!
781 : !! OUTPUT
782 : !! or2vals(skw%nr)=||R||**2
783 : !!
784 : !! SOURCE
785 :
786 6 : subroutine find_rstar_gen(skw, cryst, nrwant, rmax, or2vals, comm)
787 :
788 : !Arguments ------------------------------------
789 : !scalars
790 : type(skw_t),intent(inout) :: skw
791 : type(crystal_t),intent(in) :: cryst
792 : integer,intent(in) :: nrwant,comm
793 : !arrays
794 : integer,intent(in) :: rmax(3)
795 : real(dp),allocatable,intent(out) :: or2vals(:)
796 :
797 : !Local variables-------------------------------
798 : !scalars
799 : integer :: cnt,nstars,i1,i2,i3,msize,ir,nsh,ish,ss,ee,nst,ierr,nprocs,my_rank,ii
800 : real(dp) :: r2_prev
801 : !character(len=500) :: msg
802 : !arrays
803 6 : integer,allocatable :: iperm(:),rtmp(:,:),rgen(:,:),r2sh(:),shlim(:),sh_start(:),sh_stop(:)
804 6 : integer,allocatable :: recvcounts(:),displs(:),recvbuf(:,:)
805 6 : real(dp),allocatable :: r2tmp(:),cnorm(:)
806 :
807 : ! *********************************************************************
808 :
809 6 : nprocs = xmpi_comm_size(comm); my_rank = xmpi_comm_rank(comm)
810 :
811 24 : msize = product(2*rmax + 1)
812 18 : ABI_MALLOC(rtmp, (3, msize))
813 18 : ABI_MALLOC(r2tmp, (msize))
814 :
815 6 : cnt = 0
816 114 : do i3=-rmax(3),rmax(3)
817 2120 : do i2=-rmax(2),rmax(2)
818 40766 : do i1=-rmax(1),rmax(1)
819 38652 : cnt = cnt + 1
820 154608 : rtmp(:, cnt) = [i1,i2,i3]
821 775046 : r2tmp(cnt) = dot_product(rtmp(:,cnt), matmul(cryst%rmet, rtmp(:,cnt)))
822 : end do
823 : end do
824 : end do
825 :
826 : ! Sort r2tmp
827 18 : ABI_MALLOC(iperm, (msize))
828 115974 : iperm = [(i1, i1=1,msize)]
829 6 : call sort_dp(msize, r2tmp, iperm, tol12)
830 :
831 : ! Find R-points generating the stars.
832 12 : ABI_MALLOC(rgen, (3, msize))
833 38658 : do ir=1,msize
834 154614 : rgen(:,ir) = rtmp(:,iperm(ir))
835 : end do
836 154620 : rtmp = rgen
837 6 : ABI_FREE(iperm)
838 :
839 12 : ABI_MALLOC(r2sh, (msize)) ! Correspondence between R and the shell index.
840 18 : ABI_MALLOC(shlim, (msize+1)) ! For each shell, the index of the initial G-vector.
841 6 : nsh = 1; r2sh(1) = 1; shlim(1) = 1; r2_prev = zero
842 38652 : do ir=2,msize
843 38646 : if (abs(r2tmp(ir) - r2_prev) > r2tmp(ir) * tol8) then
844 1267 : r2_prev = r2tmp(ir); nsh = nsh + 1; shlim(nsh) = ir
845 : !write(std_out,*)"nsh: ",shlim(nsh) - shlim(nsh-1)
846 : end if
847 38652 : r2sh(ir) = nsh
848 : end do
849 6 : shlim(nsh+1) = msize + 1
850 6 : ABI_FREE(r2tmp)
851 6 : ABI_FREE(r2sh)
852 :
853 : !call get_irredg(msize, skw%ptg_nsym, +1, cryst%rprimd, skw%ptg_symrel, rtmp, nstars, rgen, cnorm)
854 : !write(66,*)nstars; do ish=1,nstars; write(66,*)rgen(:,ish); end do
855 :
856 : ! Distribute shells among processor so that we can parallelize the search algorithm.
857 : ! Each proc works on a contigous block of shells, then we have to gather the results.
858 18 : ABI_MALLOC(sh_start, (0:nprocs-1))
859 12 : ABI_MALLOC(sh_stop, (0:nprocs-1))
860 6 : call xmpi_split_work2_i4b(nsh, nprocs, sh_start, sh_stop)
861 :
862 12 : ABI_MALLOC(cnorm, (msize))
863 6 : nstars = 0
864 1279 : do ish=sh_start(my_rank),sh_stop(my_rank)
865 1273 : ss = shlim(ish); ee = shlim(ish+1) - 1; msize = ee - ss + 1
866 : call get_irredg(msize, skw%ptg_nsym, + 1, cryst%rprimd, skw%ptg_symrel, rtmp(:,ss:), &
867 1273 : nst, rgen(:,nstars+1:), cnorm(nstars+1:))
868 1279 : nstars = nstars + nst
869 : end do
870 :
871 6 : ABI_FREE(cnorm)
872 6 : ABI_FREE(sh_start)
873 6 : ABI_FREE(sh_stop)
874 6 : ABI_FREE(rtmp)
875 6 : ABI_FREE(shlim)
876 :
877 6 : if (nprocs > 1) then
878 : ! Collect star functions.
879 0 : ABI_MALLOC(recvcounts, (nprocs))
880 0 : recvcounts = 0; recvcounts(my_rank+1) = 3 * nstars
881 0 : call xmpi_sum(recvcounts, comm, ierr)
882 0 : ABI_MALLOC(displs, (nprocs))
883 0 : displs(1) = 0
884 0 : do ii=2,nprocs
885 0 : displs(ii) = sum(recvcounts(:ii-1))
886 : end do
887 0 : call xmpi_sum(nstars, nst, comm, ierr) ! Now nst is the total number of star functions.
888 0 : ABI_MALLOC(recvbuf, (3, nst))
889 0 : call xmpi_allgatherv(rgen, 3*nstars, recvbuf, recvcounts, displs, comm, ierr)
890 0 : ABI_FREE(recvcounts)
891 0 : ABI_FREE(displs)
892 0 : nstars = nst
893 0 : rgen(:,1:nstars) = recvbuf
894 0 : ABI_FREE(recvbuf)
895 : end if
896 : !if (my_rank == 0) then
897 : ! write(67,*)"nstars",nstars,"nsh",nsh; do ish=1,nstars; write(67,*)rgen(:,ish); end do
898 : !end if
899 :
900 : ! Store rpts and compute ||R||**2.
901 6 : skw%nr = min(nstars, nrwant)
902 6 : if (allocated(skw%rpts)) then
903 0 : ABI_FREE(skw%rpts)
904 : end if
905 18 : ABI_MALLOC(skw%rpts, (3, skw%nr))
906 1352 : skw%rpts = rgen(:,1:skw%nr)
907 18 : ABI_MALLOC(or2vals, (skw%nr))
908 341 : do ir=1,skw%nr
909 6706 : or2vals(ir) = dot_product(skw%rpts(:,ir), matmul(cryst%rmet, skw%rpts(:,ir)))
910 : end do
911 :
912 6 : ABI_FREE(rgen)
913 :
914 6 : end subroutine find_rstar_gen
915 : !!***
916 :
917 94907 : end module m_skw
918 : !!***
|
