Branch data 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 : 14 : 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 : 7 : integer,allocatable :: ipiv(:)
188 : : real(dp) :: list2(2)
189 : 7 : real(dp),allocatable :: r2vals(:),inv_rhor(:),oeig(:)
190 : 7 : complex(dp),allocatable :: srk(:,:),hmat(:,:),lambda(:,:,:),work(:)
191 : : ! *********************************************************************
192 : :
193 [ - + ]: 7 : ABI_CHECK(nkpt > 1, sjoin("nkpt must be > 1 but got:", itoa(nkpt)))
194 : :
195 : 7 : call cwtime(cpu_tot, wall_tot, gflops_tot, "start")
196 : :
197 : 7 : nprocs = xmpi_comm_size(comm); my_rank = xmpi_comm_rank(comm)
198 : :
199 : : ! Get slice of bands to be treated.
200 [ + + + - : 21 : new%band_block = band_block; if (all(band_block == 0)) new%band_block = [1, nband]
- + - + -
- ]
201 : 7 : bstart = new%band_block(1); bstop = new%band_block(2); bcount = bstop - bstart + 1
202 : 7 : new%cplex = cplex; new%nkpt = nkpt; new%nsppol = nsppol; new%bcount = bcount
203 : :
204 : : ! Get point group operations.
205 [ - + - + ]: 7 : 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 : 7 : lpratio = int(abs(params(1)))
211 [ - + ]: 7 : ABI_CHECK(lpratio > 0, "lpratio must be > 0")
212 [ + + ]: 28 : rmax = nint((one + (lpratio * new%nkpt * new%ptg_nsym) / two) ** third)
213 [ + - ]: 7 : if (new%has_inversion) then
214 [ + + ]: 28 : rmax = nint((one + (lpratio * new%nkpt * new%ptg_nsym / 2) / two) ** third)
215 : : end if
216 : 7 : nrwant = lpratio * new%nkpt
217 : :
218 : 7 : call cwtime(cpu, wall, gflops, "start")
219 : 0 : do
220 [ - + ]: 7 : call find_rstar_gen(new, cryst, nrwant, rmax, r2vals, comm)
221 [ - + ]: 7 : 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 : 7 : nr = new%nr
231 : 7 : call cwtime_report(" find_rstar_gen", cpu, wall, gflops)
232 : :
233 [ + - ]: 7 : if (my_rank == master) call new%print(std_out)
234 : :
235 : : ! Compute (inverse) roughness function.
236 : 7 : r2min = r2vals(2)
237 [ - + + - : 21 : ABI_MALLOC(inv_rhor, (nr))
+ - + - -
+ - + ]
238 [ + + ]: 382 : do ir=1,nr
239 : 375 : r2 = r2vals(ir)
240 : 382 : 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 [ - + + - : 28 : ABI_MALLOC(srk, (nr, nkpt))
- + + - +
- + - - +
- + ]
247 [ + + ]: 82 : do ik=1,nkpt
248 : 82 : call mkstar(new, kpts(:,ik), srk(:,ik))
249 : : end do
250 : :
251 : : ! Build H(k,k') matrix (Hermitian)
252 [ + - - + : 1150 : ABI_CALLOC(hmat, (nkpt-1, nkpt-1))
+ - - + +
- + - + -
- + - + +
+ + + ]
253 : : cnt = 0
254 [ + + ]: 75 : do jj=1,nkpt-1
255 [ + + ]: 636 : do ii=1,jj
256 : : !do ii=1,nkpt-1
257 [ + - ]: 561 : cnt = cnt + 1; if (mod(cnt, nprocs) /= my_rank) cycle ! mpi parallelism.
258 [ + + ]: 64988 : do ir=2,nr
259 : : hmat(ii, jj) = hmat(ii, jj) + &
260 : 64920 : (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 : 7 : call xmpi_sum(hmat, comm, ierr)
265 : :
266 [ - + + - : 35 : ABI_MALLOC(lambda, (nkpt-1, bcount, nsppol))
- + + - -
+ + - + -
+ - - + -
+ ]
267 [ + + ]: 14 : do spin=1,nsppol
268 [ + + ]: 98 : do ib=1,bcount
269 : 84 : band = ib + bstart - 1
270 [ + + ]: 819 : 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 : 7 : call wrtout(std_out, " Solving system of linear equations to get lambda coeffients (eq. 10 of PRB 38 2721)...", do_flush=.True.)
276 : 7 : call cwtime(cpu, wall, gflops, "start")
277 [ - + + - : 21 : 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 [ + + ]: 75 : do ii=1,nkpt-1
286 : 75 : hmat(ii, ii) = real(hmat(ii, ii))
287 : : end do
288 : 7 : lwork = -1
289 [ - + ]: 7 : ABI_MALLOC(work, (1))
290 : 7 : call zhesv("U", nkpt-1, bcount*nsppol, hmat, nkpt-1, ipiv, lambda, nkpt-1, work, lwork, ierr)
291 : 7 : lwork = nint(real(work(1)))
292 : 7 : ABI_FREE(work)
293 [ - + + - : 21 : ABI_MALLOC(work, (lwork))
+ - + - -
+ - + ]
294 : 7 : call zhesv("U", nkpt-1, bcount*nsppol, hmat, nkpt-1, ipiv, lambda, nkpt-1, work, lwork, ierr)
295 [ - + ]: 7 : ABI_CHECK(ierr == 0, sjoin("ZHESV returned:", itoa(ierr)))
296 [ - + ]: 7 : ABI_FREE(work)
297 : : end if
298 : 7 : call cwtime_report(" ZHESV", cpu, wall, gflops)
299 : :
300 : : ! Compute coefficients
301 [ - + + - : 35 : ABI_MALLOC(new%coefs, (nr,bcount,nsppol))
- + + - -
+ + - + -
+ - - + -
+ - + ]
302 : :
303 [ + + ]: 14 : do spin=1,nsppol
304 [ + + ]: 98 : do ib=1,bcount
305 : 84 : band = ib + bstart - 1
306 [ + + ]: 4060 : do ir=2,nr
307 [ + + ]: 55272 : 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 [ + + ]: 4067 : 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 [ - + ]: 7 : 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 [ + + ]: 28 : new%cached_kpt = huge(one)
327 [ - + + - : 21 : ABI_MALLOC(new%cached_srk, (new%nr))
+ - + - -
+ - + -
+ ]
328 [ + + ]: 28 : new%cached_kpt_dk1 = huge(one)
329 [ - + + - : 21 : ABI_MALLOC(new%cached_srk_dk1, (new%nr, 3))
+ - + - -
+ - + -
+ ]
330 [ + + ]: 28 : new%cached_kpt_dk2 = huge(one)
331 [ - + + - : 28 : ABI_MALLOC(new%cached_srk_dk2, (new%nr, 3, 3))
+ - + - +
- - + - +
- + ]
332 : :
333 : 7 : ABI_FREE(r2vals)
334 [ - + ]: 7 : ABI_FREE(srk)
335 [ - + ]: 7 : ABI_FREE(inv_rhor)
336 [ - + ]: 7 : ABI_FREE(hmat)
337 [ - + ]: 7 : ABI_FREE(lambda)
338 [ - + ]: 7 : ABI_FREE(ipiv)
339 : :
340 : : ! Compare ab-initio data with interpolated results.
341 [ - + + - : 21 : ABI_MALLOC(oeig, (bcount))
+ - + - -
+ - + ]
342 [ + - ]: 7 : fmt = sjoin("(a,", itoa(bcount), "(es12.4))")
343 : 7 : bstop = bstart + bcount - 1
344 : 7 : mare = zero; mae_meV = zero; cnt = 0
345 : 7 : call wrtout(std_out, ch10//" Comparing ab-initio energies with SKW interpolated results...")
346 [ + + ]: 14 : do spin=1,nsppol
347 [ + + ]: 89 : do ik=1,nkpt
348 [ - + ]: 75 : cnt = cnt + 1; if (mod(cnt, nprocs) /= my_rank) cycle ! mpi parallelism.
349 : :
350 [ + + ]: 887 : do ib=1,bcount
351 : 812 : band = ib + new%band_block(1) - 1
352 : 812 : call new%eval_bks(band, kpts(:,ik), spin, oeig(ib))
353 : :
354 : 812 : adiff_meV = abs(eig(band,ik,spin) - oeig(ib)); rel_err = zero
355 [ + - ]: 812 : if (abs(eig(band,ik,spin)) > tol16) rel_err = adiff_meV / abs(eig(band,ik,spin))
356 : 812 : rel_err = 100 * rel_err; adiff_meV = adiff_meV * Ha_meV
357 : 887 : mae_meV = mae_meV + adiff_meV; mare = mare + rel_err
358 : : end do
359 : :
360 : 7 : if (prtvol > 0) then
361 [ + - + + ]: 887 : ib = imax_loc(eig(bstart:bstop,ik,spin) - oeig)
362 : 75 : rval = (eig(bstart+ib-1,ik,spin) - oeig(ib)) * Ha_meV
363 : : write(std_out,"(a,es12.4,2a)") &
364 : 75 : " SKW maxerr: ", rval, &
365 : 150 : " (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 [ - + ]: 7 : ABI_FREE(oeig)
373 : :
374 : : ! Issue warning if error too large.
375 [ + + ]: 21 : list2 = [mare, mae_meV]; call xmpi_sum(list2, comm, ierr); mare = list2(1); mae_meV = list2(2)
376 : 7 : cnt = bcount * nkpt * nsppol; mare = mare / cnt; mae_meV = mae_meV / cnt
377 : 7 : write(std_out,"(2(a,es12.4),a,/)")" MARE: ",mare, ", MAE: ", mae_meV, " (meV)"
378 [ - + ]: 7 : 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 : 7 : call cwtime_report(" skw_new", cpu_tot, wall_tot, gflops_tot, end_str=ch10)
386 : :
387 [ - + ]: 56 : 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 : 7 : 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 : 7 : write(unt,"(a)")" === Shankland-Koelling-Wood Fourier interpolation scheme ==="
417 : 7 : write(unt,"(a)")sjoin(" nsppol", itoa(skw%nsppol), ", cplex:", itoa(skw%cplex))
418 : 7 : write(unt,"(a)")sjoin(" Number of ab-initio k-points:", itoa(skw%nkpt))
419 : 7 : write(unt,"(a)")sjoin(" Number of star functions:", itoa(skw%nr))
420 : 7 : write(unt,"(a)")sjoin(" Stars/Nk ratio:", ftoa(skw%nr * one / skw%nkpt))
421 : 7 : write(unt,"(a)")sjoin(" Has spatial inversion:", yesno(skw%has_inversion))
422 : :
423 : 7 : 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 : 13594 : 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 : 13594 : ib = band - skw%band_block(1) + 1
535 [ + - - + ]: 13594 : 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 [ + + + + : 51165 : if (any(kpt /= skw%cached_kpt)) then
+ + ]
539 : 1118 : call mkstar(skw, kpt, skw%cached_srk)
540 [ + + ]: 4472 : skw%cached_kpt = kpt
541 : : end if
542 : :
543 [ + + ]: 860874 : oeig = dot_product(conjg(skw%coefs(:,ib,spin)), skw%cached_srk)
544 : :
545 : : ! TODO: Test Derivatives
546 [ - + ]: 13594 : 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 [ - + ]: 13594 : 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 : 13594 : 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 : 7 : subroutine skw_free(skw)
589 : :
590 : : !Arguments ------------------------------------
591 : : !scalars
592 : : class(skw_t),intent(inout) :: skw
593 : :
594 : : ! *********************************************************************
595 : :
596 [ + - ]: 7 : ABI_SFREE(skw%rpts)
597 [ + - ]: 7 : ABI_SFREE(skw%ptg_symrel)
598 [ + - ]: 7 : ABI_SFREE(skw%ptg_symrec)
599 [ + - ]: 7 : ABI_SFREE(skw%coefs)
600 : :
601 [ + - ]: 7 : ABI_SFREE(skw%cached_srk)
602 [ + + ]: 28 : skw%cached_kpt = huge(one)
603 [ + - ]: 7 : ABI_SFREE(skw%cached_srk_dk1)
604 [ + + ]: 28 : skw%cached_kpt_dk1 = huge(one)
605 [ + - ]: 7 : ABI_SFREE(skw%cached_srk_dk2)
606 [ + + ]: 28 : skw%cached_kpt_dk2 = huge(one)
607 : :
608 : 7 : 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 : 1193 : 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 [ + + ]: 92893 : srk = zero
645 [ + + ]: 58457 : do isym=1,skw%ptg_nsym
646 [ + + + + : 1660656 : sk = two_pi * matmul(transpose(skw%ptg_symrel(:,:,isym)), kpt)
+ + + + +
- + + ]
647 [ + + ]: 4460057 : do ir=1,skw%nr
648 [ + + ]: 17663664 : srk(ir) = srk(ir) + exp(j_dpc * dot_product(sk, skw%rpts(:,ir)))
649 : : end do
650 : : end do
651 [ + + ]: 92893 : srk = srk / skw%ptg_nsym
652 : :
653 : 1193 : 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 : 7 : 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 : 7 : integer,allocatable :: iperm(:),rtmp(:,:),rgen(:,:),r2sh(:),shlim(:),sh_start(:),sh_stop(:)
804 : 7 : integer,allocatable :: recvcounts(:),displs(:),recvbuf(:,:)
805 : 7 : real(dp),allocatable :: r2tmp(:),cnorm(:)
806 : :
807 : : ! *********************************************************************
808 : :
809 : 7 : nprocs = xmpi_comm_size(comm); my_rank = xmpi_comm_rank(comm)
810 : :
811 [ + + ]: 28 : msize = product(2*rmax + 1)
812 [ - + + - : 21 : ABI_MALLOC(rtmp, (3, msize))
+ - + - -
+ - + ]
813 [ - + + - : 21 : ABI_MALLOC(r2tmp, (msize))
+ - + - -
+ - + ]
814 : :
815 : 7 : cnt = 0
816 [ + + ]: 132 : do i3=-rmax(3),rmax(3)
817 [ + + ]: 2427 : do i2=-rmax(2),rmax(2)
818 [ + + ]: 45985 : do i1=-rmax(1),rmax(1)
819 : 43565 : cnt = cnt + 1
820 [ + + ]: 174260 : rtmp(:, cnt) = [i1,i2,i3]
821 [ + + + + : 873595 : 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 [ - + + - : 21 : ABI_MALLOC(iperm, (msize))
+ - + - -
+ - + ]
828 [ + - + + : 130716 : iperm = [(i1, i1=1,msize)]
+ - + + +
- - + - -
+ + ]
829 : 7 : call sort_dp(msize, r2tmp, iperm, tol12)
830 : :
831 : : ! Find R-points generating the stars.
832 [ - + + - : 14 : ABI_MALLOC(rgen, (3, msize))
- + - + ]
833 [ + + ]: 43572 : do ir=1,msize
834 [ + + ]: 174267 : rgen(:,ir) = rtmp(:,iperm(ir))
835 : : end do
836 [ - + + + : 174274 : rtmp = rgen
+ + ]
837 [ - + ]: 7 : ABI_FREE(iperm)
838 : :
839 [ - + + - : 14 : ABI_MALLOC(r2sh, (msize)) ! Correspondence between R and the shell index.
- + - + ]
840 [ + - - + : 21 : ABI_MALLOC(shlim, (msize+1)) ! For each shell, the index of the initial G-vector.
+ - + - +
- - + -
+ ]
841 : 7 : nsh = 1; r2sh(1) = 1; shlim(1) = 1; r2_prev = zero
842 [ + + ]: 43565 : do ir=2,msize
843 [ + + ]: 43558 : if (abs(r2tmp(ir) - r2_prev) > r2tmp(ir) * tol8) then
844 : 1447 : r2_prev = r2tmp(ir); nsh = nsh + 1; shlim(nsh) = ir
845 : : !write(std_out,*)"nsh: ",shlim(nsh) - shlim(nsh-1)
846 : : end if
847 : 43565 : r2sh(ir) = nsh
848 : : end do
849 : 7 : shlim(nsh+1) = msize + 1
850 [ - + ]: 7 : ABI_FREE(r2tmp)
851 : 7 : 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 [ - + + - : 21 : ABI_MALLOC(sh_start, (0:nprocs-1))
+ - + - -
+ - + ]
859 [ - + + - : 14 : ABI_MALLOC(sh_stop, (0:nprocs-1))
- + - + ]
860 : 7 : call xmpi_split_work2_i4b(nsh, nprocs, sh_start, sh_stop)
861 : :
862 [ - + + - : 14 : ABI_MALLOC(cnorm, (msize))
- + - + ]
863 : 7 : nstars = 0
864 [ + + ]: 1461 : do ish=sh_start(my_rank),sh_stop(my_rank)
865 : 1454 : 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 : 1454 : nst, rgen(:,nstars+1:), cnorm(nstars+1:))
868 : 1461 : nstars = nstars + nst
869 : : end do
870 : :
871 [ - + ]: 7 : ABI_FREE(cnorm)
872 : 7 : ABI_FREE(sh_start)
873 : 7 : ABI_FREE(sh_stop)
874 [ - + ]: 7 : ABI_FREE(rtmp)
875 : 7 : ABI_FREE(shlim)
876 : :
877 [ - + ]: 7 : 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 : 7 : skw%nr = min(nstars, nrwant)
902 [ - + ]: 7 : if (allocated(skw%rpts)) then
903 : 0 : ABI_FREE(skw%rpts)
904 : : end if
905 [ - + + - : 21 : ABI_MALLOC(skw%rpts, (3, skw%nr))
+ - + - -
+ - + -
+ ]
906 [ - + + + : 1514 : skw%rpts = rgen(:,1:skw%nr)
+ + ]
907 [ - + + - : 21 : ABI_MALLOC(or2vals, (skw%nr))
+ - + - -
+ - + -
+ ]
908 [ + + ]: 382 : do ir=1,skw%nr
909 [ + + + + : 7507 : or2vals(ir) = dot_product(skw%rpts(:,ir), matmul(cryst%rmet, skw%rpts(:,ir)))
+ + + - +
+ ]
910 : : end do
911 : :
912 [ - + ]: 7 : ABI_FREE(rgen)
913 : :
914 [ - + - + : 7 : end subroutine find_rstar_gen
- + ]
915 : : !!***
916 : :
917 [ - - # # : 101204 : end module m_skw
# # # # #
# # # # #
# # # # #
# # # # #
# # ]
918 : : !!***
|
