Line data Source code
1 : !!****m* ABINIT/m_kg
2 : !! NAME
3 : !! m_kg
4 : !!
5 : !! FUNCTION
6 : !! Low-level functions to operate of G-vectors.
7 : !!
8 : !! COPYRIGHT
9 : !! Copyright (C) 2008-2026 ABINIT group (DCA, XG, GMR, MT, DRH, AR)
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_kg
23 :
24 : use defs_basis
25 : use m_errors
26 : use m_abicore
27 : use m_errors
28 : use m_xmpi
29 : use m_dtset
30 :
31 : use defs_abitypes, only : MPI_type
32 : use m_fftcore, only : kpgsph, bound
33 : use m_mpinfo, only : proc_distrb_cycle
34 : use m_time, only : timab
35 :
36 : implicit none
37 :
38 : private
39 : !!***
40 :
41 : public :: getcut ! Compute cutoff G^2
42 : public :: getmpw ! Compute recommended npw from ecut, ucvol and gmet
43 : public :: mkkin ! Compute elements of kinetic energy operator in reciprocal space at a given k point
44 : public :: kpgio ! Do initialization of kg data.
45 : public :: ph1d3d ! Compute the three-dimensional phase factor $e^{i 2 \pi (k+G) cdot xred}$
46 : public :: getph ! Compute three factors of one-dimensional structure factor phase
47 : public :: kpgstr ! Derivative of kinetic energy operator in reciprocal space.
48 : public :: mkkpg ! Compute all (k+G) vectors (dp, in reduced coordinates) for given k point
49 : public :: mkpwind_k ! Make plane wave index at k point for basis at second k point
50 : public :: mkkpgcart ! Compute all (k+G) vectors (dp, in cartesian coordinates) for given k point
51 : public :: mkkin_metdqdq ! Compute the second q-gradient of the derivative of the kinetic energy operator w.r.t a metric
52 :
53 : contains
54 : !!***
55 :
56 : !!****f* m_kg/getcut
57 : !! NAME
58 : !! getcut
59 : !!
60 : !! FUNCTION
61 : !! For input kpt, fft box dim ngfft(1:3), recip space metric gmet, and kinetic energy cutoff ecut, COMPUTES:
62 : !!
63 : !! if iboxcut==0:
64 : !! gsqcut: cut-off on G^2 for "large sphere" of radius double that
65 : !! of the basis sphere corresponding to ecut
66 : !! boxcut: where boxcut == gcut(box)/gcut(sphere).
67 : !! boxcut >=2 for no aliasing.
68 : !! boxcut < 1 is wrong and halts subroutine.
69 : !! if iboxcut==1:
70 : !! gsqcut: cut-off on G^2 for "large sphere" containing the whole fft box
71 : !! boxcut: no meaning (zero)
72 : !!
73 : !! INPUTS
74 : !! ecut=kinetic energy cutoff for planewave sphere (hartree)
75 : !! gmet(3,3)=reciprocal space metric (bohr^-2)
76 : !! iboxcut=0: compute gsqcut and boxcut with boxcut>=1
77 : !! 1: compute gsqcut for boxcut=1 (sphere_cutoff=box_cutoff)
78 : !! iout=unit number for output file
79 : !! kpt(3)=input k vector (reduced coordinates--in terms of reciprocal lattice primitive translations)
80 : !! ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft
81 : !!
82 : !! OUTPUT
83 : !! boxcut=defined above (dimensionless), ratio of basis sphere
84 : !! diameter to fft box length (smallest value)
85 : !! gsqcut=Fourier cutoff on G^2 for "large sphere" of radius double
86 : !! that of the basis sphere--appropriate for charge density rho(G),
87 : !! Hartree potential, and pseudopotentials
88 : !!
89 : !! NOTES
90 : !! 2*gcut arises from rho(g)=sum g prime (psi(g primt)*psi(g prime+g))
91 : !! where psi(g) is only nonzero for |g| <= gcut).
92 : !! ecut (currently in hartree) is proportional to gcut(sphere)**2.
93 : !!
94 : !! SOURCE
95 :
96 38564 : subroutine getcut(boxcut, ecut, gmet, gsqcut, iboxcut, iout, kpt, ngfft)
97 :
98 : !Arguments ------------------------------------
99 : !scalars
100 : integer,intent(in) :: iboxcut,iout
101 : real(dp),intent(in) :: ecut
102 : real(dp),intent(out) :: boxcut,gsqcut
103 : !arrays
104 : integer,intent(in) :: ngfft(18)
105 : real(dp),intent(in) :: gmet(3,3),kpt(3)
106 :
107 : !Local variables-------------------------------
108 : !scalars
109 : integer :: plane
110 : real(dp) :: boxsq,cutrad,ecut_pw,effcut,largesq,sphsq
111 : character(len=1000) :: msg
112 : !arrays
113 : integer :: gbound(3)
114 : ! *************************************************************************
115 :
116 : ! This is to treat the case in which ecut has not been initialized e.g. for wavelet computations.
117 : ! The default for ecut is -1.0 , allowed only for wavelets calculations
118 19282 : ecut_pw=ecut
119 19282 : if(ecut<-tol8)ecut_pw=ten
120 :
121 : !gcut(box)**2=boxsq; gcut(sphere)**2=sphsq
122 : !get min. d**2 to boundary of fft box:
123 : !(gmet sets dimensions: bohr**-2)
124 : !ecut(sphere)=0.5*(2 pi)**2 * sphsq:
125 19282 : call bound(largesq,boxsq,gbound,gmet,kpt,ngfft,plane)
126 19282 : effcut=0.5_dp * (two_pi)**2 * boxsq
127 19282 : sphsq=2._dp*ecut_pw/two_pi**2
128 :
129 19282 : if (iboxcut/=0) then
130 3 : boxcut=10._dp
131 3 : gsqcut=(largesq/sphsq)*(2.0_dp*ecut)/two_pi**2
132 :
133 3 : write(msg, '(a,a,3f8.4,a,3i4,a,a,f11.3,a,a)' ) ch10,&
134 3 : ' getcut: wavevector=',kpt,' ngfft=',ngfft(1:3),ch10,&
135 6 : ' ecut(hartree)=',ecut_pw+tol8,ch10,'=> whole FFT box selected'
136 3 : if(iout/=std_out) call wrtout(iout,msg)
137 3 : call wrtout(std_out,msg)
138 : else
139 :
140 : ! Get G^2 cutoff for sphere of double radius of basis sphere
141 : ! for selecting G s for rho(G), V_Hartree(G), and V_psp(G)--
142 : ! cut off at fft box boundary or double basis sphere radius, whichever
143 : ! is smaller. If boxcut were 2, then relation would be
144 : ! $ecut_eff = (1/2) * (2 Pi Gsmall)^2 and gsqcut=4*Gsmall^2$.
145 19279 : boxcut = sqrt(boxsq/sphsq)
146 19279 : cutrad = min(2.0_dp,boxcut)
147 19279 : gsqcut = (cutrad**2)*(2.0_dp*ecut_pw)/two_pi**2
148 :
149 19279 : if(ecut>-tol8)then
150 :
151 19279 : write(msg, '(a,a,3f8.4,a,3i4,a,a,f11.3,3x,a,f10.5)' ) ch10,&
152 19279 : ' getcut: wavevector=',kpt,' ngfft=',ngfft(1:3),ch10,&
153 38558 : ' ecut(hartree)=',ecut+tol8,'=> boxcut(ratio)=',boxcut+tol8
154 19279 : if(iout/=std_out) call wrtout(iout,msg)
155 19279 : call wrtout(std_out,msg)
156 :
157 19279 : if (boxcut<1.0_dp) then
158 : write(msg, '(9a,f12.6,6a)' )&
159 0 : 'Choice of acell, ngfft, and ecut',ch10,&
160 0 : '===> basis sphere extends BEYOND FFT box !',ch10,&
161 0 : 'Recall that boxcut=Gcut(box)/Gcut(sphere) must be > 1.',ch10,&
162 0 : 'Action: try larger ngfft or smaller ecut.',ch10,&
163 0 : 'Note that ecut=effcut/boxcut**2 and effcut=',effcut+tol8,ch10,&
164 0 : 'This situation might happen when optimizing the cell parameters.',ch10,&
165 0 : 'Your starting geometry might be crazy.',ch10,&
166 0 : 'See https://wiki.abinit.org/doku.php?id=howto:troubleshooting#incorrect_initial_geometry .'
167 0 : if(iout/=std_out) call wrtout(iout,msg)
168 0 : ABI_ERROR(msg)
169 : end if
170 :
171 19279 : if (boxcut>2.2_dp) then
172 7145 : write(msg, '(a,a,a,a,a,a,a,a,a,a,a,f12.6,a,a)' ) ch10,&
173 7145 : ' getcut : COMMENT -',ch10,&
174 7145 : ' Note that boxcut > 2.2 ; recall that',' boxcut=Gcut(box)/Gcut(sphere) = 2',ch10,&
175 7145 : ' is sufficient for exact treatment of convolution.',ch10,&
176 7145 : ' Such a large boxcut is a waste : you could raise ecut',ch10,&
177 14290 : ' e.g. ecut=',effcut*0.25_dp+tol8,' Hartrees makes boxcut=2',ch10
178 7145 : if(iout/=std_out) call wrtout(iout,msg)
179 7145 : call wrtout(std_out,msg)
180 : end if
181 :
182 19279 : if (boxcut<1.5_dp) then
183 41 : write(msg, '(15a)' ) ch10,&
184 41 : ' getcut : WARNING -',ch10,&
185 41 : ' Note that boxcut < 1.5; this usually means',ch10,&
186 41 : ' that the forces are being fairly strongly affected by',' the smallness of the FFT box.',ch10,&
187 41 : ' Be sure to test with larger ngfft(1:3) values.',ch10,&
188 41 : ' This situation might happen when optimizing the cell parameters.',ch10,&
189 41 : ' Your starting geometry might be crazy.',ch10,&
190 82 : ' See https://wiki.abinit.org/doku.php?id=howto:troubleshooting#incorrect_initial_geometry .'
191 41 : if(iout/=std_out) call wrtout(iout,msg)
192 41 : call wrtout(std_out,msg)
193 : end if
194 :
195 : end if
196 :
197 : end if ! iboxcut
198 :
199 19282 : end subroutine getcut
200 : !!***
201 :
202 : !!****f* m_kg/getmpw
203 : !! NAME
204 : !! getmpw
205 : !!
206 : !! FUNCTION
207 : !! From input ecut, combined with ucvol and gmet, compute recommended mpw
208 : !! mpw is the maximum number of plane-waves in the wave-function basis
209 : !! for one processor of the WF group
210 : !!
211 : !! INPUTS
212 : !! ecut=plane-wave cutoff energy in Hartrees
213 : !! exchn2n3d=if n1, n2 and n3 are exchanged
214 : !! gmet(3,3)=reciprocal space metric (bohr**-2).
215 : !! istwfk(nkpt)=input option parameter that describes the storage of wfs
216 : !! kptns(3,nkpt)=real(dp) array for k points (normalisation is already taken into account)
217 : !! mpi_enreg=information about MPI parallelization
218 : !! nkpt=integer number of k points in the calculation
219 : !!
220 : !! OUTPUT
221 : !! mpw=maximal number of plane waves over all k points of the processor
222 : !! (for one processor of the WF group)
223 : !!
224 : !! SOURCE
225 :
226 12525 : subroutine getmpw(ecut,exchn2n3d,gmet,istwfk,kptns,mpi_enreg,mpw,nkpt)
227 :
228 : !Arguments ------------------------------------
229 : !scalars
230 : integer,intent(in) :: exchn2n3d,nkpt
231 : integer,intent(out) :: mpw
232 : real(dp),intent(in) :: ecut
233 : type(MPI_type),intent(inout) :: mpi_enreg
234 : !arrays
235 : integer,intent(in) :: istwfk(nkpt)
236 : real(dp),intent(in) :: gmet(3,3),kptns(3,nkpt)
237 :
238 : !Local variables-------------------------------
239 : !scalars
240 : integer :: ikpt,istwf_k,npw
241 : ! integer :: npwwrk,pad=50
242 : ! real(dp) :: scale=1.3_dp
243 : character(len=500) :: msg
244 : !arrays
245 12525 : integer,allocatable :: kg(:,:)
246 : real(dp) :: kpoint(3)
247 : ! *************************************************************************
248 :
249 : !An upper bound for mpw, might be obtained as follows
250 : !the average number of plane-waves in the cutoff sphere is
251 : !npwave = (2*ecut)**(3/2) * ucvol / (6*pi**2)
252 : !the upper bound is calculated as
253 : !npwwrk = int(scale * npwave) + pad
254 : !rescale so an upper bound
255 : !npwave=nint(ucvol*(2.0_dp*ecut)**1.5_dp/(6.0_dp*pi**2))
256 : !npwwrk=nint(dble(npwave)*scale)+pad
257 :
258 12525 : ABI_MALLOC(kg,(3,100))
259 :
260 : !set mpw to zero, as needed for only counting in kpgsph
261 12525 : mpw = 0
262 :
263 : !Might be parallelized over k points ? !
264 477100 : do ikpt = 1,nkpt
265 : ! Do computation of G sphere, returning npw
266 1858300 : kpoint(:)=kptns(:,ikpt)
267 464575 : istwf_k=istwfk(ikpt)
268 464575 : call kpgsph(ecut,exchn2n3d,gmet,0,ikpt,istwf_k,kg,kpoint,0,mpi_enreg,0,npw)
269 477100 : mpw = max(npw,mpw)
270 : end do
271 :
272 12525 : write(msg,'(a,i0)') ' getmpw: optimal value of mpw= ',mpw
273 12525 : call wrtout(std_out,msg)
274 :
275 12525 : ABI_FREE(kg)
276 :
277 12525 : end subroutine getmpw
278 : !!***
279 :
280 : !!****f* m_kg/mkkin
281 : !! NAME
282 : !! mkkin
283 : !!
284 : !! FUNCTION
285 : !! compute elements of kinetic energy operator in reciprocal space at a given k point
286 : !!
287 : !! INPUTS
288 : !! ecut=cut-off energy for plane wave basis sphere (Ha)
289 : !! ecutsm=smearing energy for plane wave kinetic energy (Ha)
290 : !! effmass_free=effective mass for electrons (1. in common case)
291 : !! gmet(3,3)=reciprocal lattice metric tensor ($\textrm{Bohr}^{-2}$)
292 : !! idir1 = 1st direction of the derivative (if 1 <= idir1 <= 3, not used otherwise)
293 : !! idir2 = 2st direction of the derivative (if 1 <= idir1,idir2 <= 3, not used otherwise))
294 : !! kg(3,npw)=integer coordinates of planewaves in basis sphere.
295 : !! kpt(3)=reduced coordinates of k point
296 : !! npw=number of plane waves at kpt.
297 : !! vecpot (optional) = vector potential used in case of RT-TDDFT with electric field
298 : !!
299 : !! OUTPUT
300 : !! kinpw(npw)=(modified) kinetic energy (or derivative) for each plane wave (Hartree)
301 : !!
302 : !! NOTES
303 : !! Usually, the kinetic energy expression is $(1/2) (2 \pi)^2 (k+G)^2 $
304 : !! However, the present implementation allows for a modification
305 : !! of this kinetic energy, in order to obtain smooth total energy
306 : !! curves with respect to the cut-off energy or the cell size and shape.
307 : !! Thus the usual expression is kept if it is lower then ecut-ecutsm,
308 : !! zero is returned beyond ecut, and in between, the kinetic
309 : !! energy is DIVIDED by a smearing factor (to make it infinite at the
310 : !! cut-off energy). The smearing factor is $x^2 (3-2x)$, where
311 : !! x = (ecut- unmodified energy)/ecutsm.
312 : !! This smearing factor is also used to derived a modified kinetic
313 : !! contribution to stress, in another routine (forstrnps.f)
314 : !! If a vector potential is given then the expression also includes
315 : !! its contributions so the kinetic energy operator is given by
316 : !! $(1/2) (2 \pi)^2 (k+G)^2 + (2 \pi) A\cdot(k+G) + (1/2) A^2$
317 : !!
318 : !! Also, in order to break slightly the symmetry between axes, that causes
319 : !! sometimes a degeneracy of eigenvalues and do not allow to obtain
320 : !! the same results on different machines, there is a modification
321 : !! by one part over 1.0e12 of the metric tensor elements (1,1) and (3,3)
322 : !!
323 : !! SOURCE
324 :
325 3651984 : subroutine mkkin(ecut, ecutsm, effmass_free, gmet, kg, kinpw, kpt, npw, idir1, idir2, vecpot)
326 :
327 : !Arguments ------------------------------------
328 : !scalars
329 : integer,intent(in) :: npw, idir1,idir2
330 : real(dp),intent(in) :: ecut,ecutsm,effmass_free
331 : !arrays
332 : integer,intent(in) :: kg(3,npw)
333 : real(dp),intent(in) :: gmet(3,3),kpt(3)
334 : real(dp),intent(out) :: kinpw(npw)
335 : real(dp),optional,intent(in) :: vecpot(3)
336 :
337 : !Local variables-------------------------------
338 : !scalars
339 : integer :: ig,order
340 : real(dp),parameter :: break_symm=1.0d-11
341 : real(dp) :: ecutsm_inv,fsm,gpk1,gpk2,gpk3,htpisq,kinetic,kpg2,dkpg2,xx
342 : real(dp) :: d1kpg2,d2kpg2,ddfsm,dfsm
343 : real(dp) :: akpg,asq
344 : logical :: l_vecpot
345 : !arrays
346 : real(dp) :: gmet_break(3,3) !, tsec(2)
347 : ! *************************************************************************
348 :
349 3651984 : l_vecpot = .false.
350 3651984 : if (present(vecpot)) then
351 27380 : if (abs(vecpot(1))>tol12 .or. abs(vecpot(2))>tol12 .or. abs(vecpot(3))>tol12) l_vecpot = .true.
352 : end if
353 :
354 : ! htpisq is (1/2) (2 Pi) **2:
355 3651984 : htpisq=0.5_dp*(two_pi)**2
356 :
357 3651984 : ecutsm_inv=0.0_dp
358 3651984 : if(ecutsm>1.0d-20)ecutsm_inv=1/ecutsm
359 :
360 3651984 : gmet_break(:,:)=gmet(:,:)
361 3651984 : gmet_break(1,1)=(1.0_dp+break_symm)*gmet(1,1)
362 3651984 : gmet_break(3,3)=(1.0_dp-break_symm)*gmet(3,3)
363 :
364 3651984 : order=0 ! Compute the kinetic operator
365 3651984 : if (idir1>0.and.idir1<4) then
366 350642 : order=1 ! Compute the 1st derivative of the kinetic operator
367 350642 : if (idir2>0.and.idir2<4) then
368 11376 : order=2 ! Compute the 2nd derivative of the kinetic operator
369 : end if
370 : end if
371 :
372 3651984 : if (l_vecpot) then
373 : ! A^2
374 : asq=( gmet_break(1,1)*vecpot(1)*vecpot(1) + &
375 : gmet_break(2,2)*vecpot(2)*vecpot(2) + &
376 : gmet_break(3,3)*vecpot(3)*vecpot(3) + &
377 : 2.0_dp*( &
378 : gmet_break(1,2)*vecpot(1)*vecpot(2) + &
379 : gmet_break(1,3)*vecpot(1)*vecpot(3) + &
380 15296 : gmet_break(2,3)*vecpot(2)*vecpot(3) ) )
381 : end if
382 :
383 : !$OMP PARALLEL DO PRIVATE(dkpg2,d1kpg2,d2kpg2,gpk1,gpk2,gpk3,ig,kinetic,kpg2,xx,fsm,dfsm,ddfsm,akpg) &
384 : !$OMP SHARED(kinpw,ecut,ecutsm,ecutsm_inv) &
385 : !$OMP SHARED(gmet_break,htpisq,idir1,idir2,kg,kpt,npw,vecpot,asq)
386 794781932 : do ig=1,npw
387 791129948 : gpk1=dble(kg(1,ig))+kpt(1)
388 791129948 : gpk2=dble(kg(2,ig))+kpt(2)
389 791129948 : gpk3=dble(kg(3,ig))+kpt(3)
390 : kpg2=htpisq*&
391 : & ( gmet_break(1,1)*gpk1**2+ &
392 : & gmet_break(2,2)*gpk2**2+ &
393 : & gmet_break(3,3)*gpk3**2 &
394 : & +2.0_dp*(gpk1*gmet_break(1,2)*gpk2+&
395 : & gpk1*gmet_break(1,3)*gpk3+ &
396 791129948 : & gpk2*gmet_break(2,3)*gpk3 ) )
397 791129948 : if (l_vecpot) then
398 : ! A.(k+G)
399 : akpg=( gmet_break(1,1)*vecpot(1)*gpk1 + &
400 : gmet_break(2,2)*vecpot(2)*gpk2 + &
401 : gmet_break(3,3)*vecpot(3)*gpk3 + &
402 : gmet_break(1,2)*(vecpot(1)*gpk2 + vecpot(2)*gpk1) + &
403 : gmet_break(1,3)*(vecpot(1)*gpk3 + vecpot(3)*gpk1) + &
404 2131560 : gmet_break(2,3)*(vecpot(2)*gpk3 + vecpot(3)*gpk2) )
405 : end if
406 743637965 : select case (order)
407 : case (0)
408 743637965 : kinetic=kpg2
409 743637965 : if (l_vecpot) kinetic=kinetic+two_pi*akpg+0.5_dp*asq
410 : case (1)
411 : dkpg2=htpisq*2.0_dp*&
412 45875003 : & (gmet_break(idir1,1)*gpk1+gmet_break(idir1,2)*gpk2+gmet_break(idir1,3)*gpk3)
413 45875003 : kinetic=dkpg2
414 45875003 : if (l_vecpot) kinetic=kinetic+two_pi*(gmet_break(idir1,1)*vecpot(1) + &
415 : & gmet_break(idir1,2)*vecpot(2) + &
416 0 : & gmet_break(idir1,3)*vecpot(3) )
417 : case (2)
418 1616980 : dkpg2=htpisq*2.0_dp*gmet_break(idir1,idir2)
419 791129948 : kinetic=dkpg2
420 : end select
421 :
422 791129948 : if(kpg2>ecut-ecutsm)then
423 8602210 : if(kpg2>ecut-tol12)then
424 1886422 : if(order==0) then
425 : ! Will filter the wavefunction, based on this value, in cgwf.f, getghc.f and precon.f
426 : kinetic=huge(zero)*1.d-10
427 : else
428 : ! The wavefunction has been filtered: no derivative
429 0 : kinetic=0
430 : end if
431 : else
432 6715788 : if(order==0) then
433 5548263 : xx=max( (ecut-kpg2)*ecutsm_inv , 1.0d-20)
434 : else
435 1167525 : xx=(ecut-kpg2)*ecutsm_inv
436 : end if
437 6715788 : if(order==2) then
438 : d1kpg2=htpisq*2.0_dp*&
439 25414 : & (gmet_break(idir1,1)*gpk1+gmet_break(idir1,2)*gpk2+gmet_break(idir1,3)*gpk3)
440 : d2kpg2=htpisq*2.0_dp*&
441 25414 : & (gmet_break(idir2,1)*gpk1+gmet_break(idir2,2)*gpk2+gmet_break(idir2,3)*gpk3)
442 : end if
443 : ! This kinetic cutoff smoothing function and its xx derivatives
444 : ! were produced with Mathematica and the fortran code has been
445 : ! numerically checked against Mathematica.
446 6715788 : fsm=1.0_dp/(xx**2*(3+xx*(1+xx*(-6+3*xx))))
447 6715788 : if(order>0) dfsm=-3.0_dp*(-1+xx)**2*xx*(2+5*xx)*fsm**2
448 1167525 : if(order>1) ddfsm=6.0_dp*xx**2*(9+xx*(8+xx*(-52+xx*(-3+xx*(137+xx*(-144+45*xx))))))*fsm**3
449 5548263 : select case (order)
450 : case (0)
451 5548263 : kinetic=kpg2*fsm
452 : case (1)
453 1142111 : kinetic=dkpg2*(fsm-ecutsm_inv*kpg2*dfsm)
454 : case (2)
455 : kinetic=dkpg2*fsm&
456 : & -2.0_dp*d1kpg2*dfsm*ecutsm_inv*d2kpg2&
457 : & +kpg2*ddfsm*(ecutsm_inv**2)*d1kpg2*d2kpg2&
458 6715788 : & -kpg2*dfsm*ecutsm_inv*dkpg2
459 : end select
460 : end if
461 : end if
462 794781932 : kinpw(ig)=kinetic/effmass_free
463 : end do
464 : !$OMP END PARALLEL DO
465 :
466 3651984 : end subroutine mkkin
467 : !!***
468 :
469 : !!****f* m_kg/kpgio
470 : !! NAME
471 : !! kpgio
472 : !!
473 : !! FUNCTION
474 : !! Do initialization of kg data.
475 : !!
476 : !! INPUTS
477 : !! ecut=kinetic energy planewave cutoff (hartree)
478 : !! exchn2n3d=if 1, n2 and n3 are exchanged
479 : !! gmet(3,3)=reciprocal space metric (bohr^-2)
480 : !! istwfk(nkpt)=input option parameter that describes the storage of wfs
481 : !! kptns(3,nkpt)=reduced coords of k points
482 : !! mkmem =number of k points treated by this node.
483 : !! character(len=4) : mode_paral=either 'COLL' or 'PERS', tells whether
484 : !! the loop over k points must be done by all processors or not,
485 : !! in case of parallel execution.
486 : !! mpi_enreg=information about MPI parallelization
487 : !! mpw=maximum number of planewaves as dimensioned in calling routine
488 : !! nband(nkpt*nsppol)=number of bands at each k point
489 : !! nkpt=number of k points
490 : !! nsppol=1 for unpolarized, 2 for polarized
491 : !!
492 : !! OUTPUT
493 : !! npwarr(nkpt)=array holding npw for each k point, taking into account
494 : !! the effect of istwfk, and the spreading over processors
495 : !! npwtot(nkpt)=array holding the total number of plane waves for each k point,
496 : !! kg(3,mpw*mkmem)=dimensionless coords of G vecs in basis sphere at k point
497 : !!
498 : !! NOTES
499 : !! Note that in case of band parallelism, the number of spin-up
500 : !! and spin-down bands must be equal at each k point.
501 : !!
502 : !! SOURCE
503 :
504 11593 : subroutine kpgio(ecut,exchn2n3d,gmet,istwfk,kg,kptns,mkmem,nband,nkpt,&
505 11593 : & mode_paral,mpi_enreg,mpw,npwarr,npwtot,nsppol)
506 :
507 : !Arguments ------------------------------------
508 : !scalars
509 : integer,intent(in) :: exchn2n3d,mkmem,mpw,nkpt,nsppol
510 : real(dp),intent(in) :: ecut
511 : character(len=4),intent(in) :: mode_paral
512 : type(MPI_type),intent(inout) :: mpi_enreg
513 : !arrays
514 : integer,intent(in) :: istwfk(nkpt),nband(nkpt*nsppol)
515 : integer,intent(out) :: kg(3,mpw*mkmem),npwarr(nkpt),npwtot(nkpt)
516 : real(dp),intent(in) :: gmet(3,3),kptns(3,nkpt)
517 :
518 : !Local variables-------------------------------
519 : !scalars
520 : integer :: ierr,ikg,ikpt,istwf_k,me,nband_down,nband_k,npw1
521 : logical :: test_npw
522 : character(len=500) :: msg
523 : !arrays
524 : real(dp) :: kpoint(3)
525 : ! *************************************************************************
526 :
527 : !Define me
528 11593 : me=mpi_enreg%me_kpt
529 :
530 11593 : if((mpi_enreg%paralbd==1) .and. (mode_paral=='PERS')) then
531 5510 : if(nsppol==2)then
532 4727 : do ikpt=1,nkpt
533 4332 : nband_k=nband(ikpt)
534 4332 : nband_down=nband(ikpt+nkpt)
535 4727 : if(nband_k/=nband_down)then
536 260 : write(msg,'(a,a,a,a,a,a,a,a,i4,a,i4,a,a,a,i4,a,a,a)')ch10,&
537 260 : ' kpgio: ERROR -',ch10,&
538 260 : ' Band parallel case, one must have same number',ch10,&
539 260 : ' of spin up and spin down bands, but input is :',ch10,&
540 260 : ' nband(up)=',nband_k,', nband(down)=',nband_down,',',ch10,&
541 260 : ' for ikpt=',ikpt,'.',ch10,&
542 520 : ' Action: correct nband in your input file.'
543 : ! MG: Tests v3(10,11,17) and v6(67) fail if this test is enabled
544 : ! call wrtout(std_out,msg,mode_paral)
545 : end if
546 : end do
547 : end if
548 : end if
549 441270 : npwarr(:)=0
550 441270 : npwtot(:)=0
551 :
552 223093905 : kg=0
553 11593 : ikg=0
554 : !Find (k+G) sphere for each k.
555 :
556 441270 : do ikpt=1,nkpt
557 :
558 429677 : nband_k = nband(ikpt)
559 :
560 429677 : if(mode_paral=='PERS')then
561 429599 : if(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,-1,me)) cycle
562 : end if
563 :
564 1405364 : kpoint(:)=kptns(:,ikpt)
565 351341 : istwf_k=istwfk(ikpt)
566 351341 : call kpgsph(ecut,exchn2n3d,gmet,ikg,ikpt,istwf_k,kg,kpoint,mkmem,mpi_enreg,mpw,npw1)
567 :
568 351341 : test_npw=.true.
569 : if (xmpi_paral==1)then
570 351341 : if (mode_paral=='PERS')then
571 2827645 : test_npw=(minval(mpi_enreg%proc_distrb(ikpt,1:nband_k,1:nsppol))==me)
572 : end if
573 : end if
574 351341 : if (test_npw) npwarr(ikpt)=npw1
575 :
576 : ! Make sure npw < nband never happens:
577 : ! if (npw1<nband(ikpt)) then
578 : ! write(msg, '(a,a,a,a,i5,a,3f8.4,a,a,i10,a,i10,a,a,a,a)' )ch10,&
579 : ! & ' kpgio : ERROR -',ch10,&
580 : ! & ' At k point number',ikpt,' k=',(kptns(ierr,ikpt),ierr=1,3),ch10,&
581 : ! & ' npw=',npw1,' < nband=',nband(ikpt),ch10,&
582 : ! & ' Indicates not enough planewaves for desired number of bands.',ch10,&
583 : ! & ' Action: change either ecut or nband in input file.'
584 : ! ABI_ERROR(msg)
585 : ! end if
586 :
587 : ! Find boundary of G sphere for efficient zero padding,
588 : ! Shift to next section of each array kg
589 792611 : ikg=ikg+npw1
590 : end do ! End of the loop over k points
591 :
592 : ! TODO: this fails on some platforms if nproc > nkpt
593 11593 : if(mode_paral == 'PERS') then
594 11515 : call xmpi_sum(npwarr,mpi_enreg%comm_kpt,ierr)
595 : end if
596 :
597 : !if (mpi_enreg%nproc>1) call wrtout(std_out,' kpgio: loop on k-points done in parallel','COLL')
598 :
599 : !XG030513 MPIWF : now, one should sum npwarr over all processors
600 : !of the WF group, to get npwtot (to be spread on all procs of the WF group
601 441270 : npwtot(:)=npwarr(:)
602 :
603 : !Taking into account istwfk
604 441270 : do ikpt=1,nkpt
605 441270 : if(istwfk(ikpt)>1)then
606 3147 : if(istwfk(ikpt)==2)then
607 1451 : npwtot(ikpt)=2*npwtot(ikpt)-1
608 : else
609 1696 : npwtot(ikpt)=2*npwtot(ikpt)
610 : end if
611 : end if
612 : end do
613 :
614 11593 : end subroutine kpgio
615 : !!***
616 :
617 : !!****f* m_kg/ph1d3d
618 : !! NAME
619 : !! ph1d3d
620 : !!
621 : !! FUNCTION
622 : !! Compute the three-dimensional phase factor $e^{i 2 \pi (k+G) cdot xred}$
623 : !! from the three one-dimensional factors, the k point coordinates,
624 : !! and the atom positions, for all planewaves which fit in the fft box.
625 : !!
626 : !! INPUTS
627 : !! iatom, jatom= bounds of atom indices in ph1d for which ph3d has to be computed
628 : !! kg_k(3,npw_k)=reduced planewave coordinates.
629 : !! matblk= dimension of ph3d
630 : !! natom= dimension of ph1d
631 : !! npw=number of plane waves
632 : !! n1,n2,n3=dimensions of fft box (ngfft(3)).
633 : !! phkxred(2,natom)=phase factors exp(2 pi k.xred)
634 : !! ph1d(2,(2*n1+1)*natom+(2*n2+1)*natom+(2*n3+1)*natom)=exp(2Pi i G xred) for
635 : !! vectors (Gx,0,0), (0,Gy,0) and (0,0,Gz)
636 : !! with components ranging from -nj <= Gj <= nj
637 : !!
638 : !! OUTPUT
639 : !! ph3d(2,npw_k,matblk)=$e^{2 i \pi (k+G) cdot xred}$ for vectors (Gx,Gy,Gz),
640 : !! and for atoms in the range iatom to jatom with respect to ph1d
641 : !!
642 : !! SOURCE
643 :
644 8172161 : subroutine ph1d3d(iatom, jatom, kg_k, matblk, natom, npw_k, n1, n2, n3, phkxred, ph1d, ph3d)
645 :
646 : !Arguments ------------------------------------
647 : !scalars
648 : integer,intent(in) :: iatom,jatom,matblk,n1,n2,n3,natom,npw_k
649 : !arrays
650 : integer,intent(in) :: kg_k(3,npw_k)
651 : real(dp),intent(in) :: ph1d(2,(2*n1+1+2*n2+1+2*n3+1)*natom)
652 : real(dp),intent(in) :: phkxred(2,natom)
653 : real(dp),intent(out) :: ph3d(2,npw_k,matblk)
654 :
655 : !Local variables-------------------------------
656 : !scalars
657 : integer :: i1,ia,iatblk,ig,kg1,kg2,kg3,shift1,shift2,shift3
658 : real(dp) :: ph12i,ph12r,ph1i,ph1r,ph2i,ph2r,ph3i,ph3r,phkxi,phkxr
659 : character(len=500) :: msg
660 : !arrays
661 8172161 : real(dp),allocatable :: ph1kxred(:,:)
662 : ! *************************************************************************
663 :
664 8172161 : if(matblk-1 < jatom-iatom)then
665 : write(msg,'(a,a,a,a,a,i0,a,a,i0,a,i0,a)')&
666 0 : 'Input natom-1 must be larger or equal to jatom-iatom,',ch10,&
667 0 : 'while their value is : ',ch10,&
668 0 : 'natom-1 = ',natom-1,ch10,&
669 0 : 'jatom=',jatom,', iatom=',iatom,'.'
670 0 : ABI_BUG(msg)
671 : end if
672 :
673 24516483 : ABI_MALLOC(ph1kxred,(2,-n1:n1))
674 :
675 : ! ia runs from iatom to jatom
676 18450831 : do ia=iatom,jatom
677 :
678 : ! iatblk runs from 1 to matblk
679 10278670 : iatblk=ia-iatom+1
680 10278670 : shift1=1+n1+(ia-1)*(2*n1+1)
681 10278670 : shift2=1+n2+(ia-1)*(2*n2+1)+natom*(2*n1+1)
682 10278670 : shift3=1+n3+(ia-1)*(2*n3+1)+natom*(2*n1+1+2*n2+1)
683 : ! Compute product of phkxred by phase for the first component of G vector
684 10278670 : phkxr=phkxred(1,ia)
685 10278670 : phkxi=phkxred(2,ia)
686 338439926 : do i1=-n1,n1
687 328161256 : ph1kxred(1,i1)=ph1d(1,i1+shift1)*phkxr-ph1d(2,i1+shift1)*phkxi
688 338439926 : ph1kxred(2,i1)=ph1d(2,i1+shift1)*phkxr+ph1d(1,i1+shift1)*phkxi
689 : end do
690 :
691 : ! Compute tri-dimensional phase factor
692 : !$OMP PARALLEL DO PRIVATE(ig,ph1r,ph1i,ph2r,ph2i,ph3r,ph3i,ph12r,ph12i,kg1,kg2,kg3)
693 2394698644 : do ig=1,npw_k
694 2376247813 : kg1=kg_k(1,ig)
695 2376247813 : kg2=kg_k(2,ig)+shift2
696 2376247813 : kg3=kg_k(3,ig)+shift3
697 2376247813 : ph1r=ph1kxred(1,kg1)
698 2376247813 : ph1i=ph1kxred(2,kg1)
699 2376247813 : ph2r=ph1d(1,kg2)
700 2376247813 : ph2i=ph1d(2,kg2)
701 2376247813 : ph3r=ph1d(1,kg3)
702 2376247813 : ph3i=ph1d(2,kg3)
703 2376247813 : ph12r=ph1r*ph2r-ph1i*ph2i
704 2376247813 : ph12i=ph1r*ph2i+ph1i*ph2r
705 2376247813 : ph3d(1,ig,iatblk)=ph12r*ph3r-ph12i*ph3i
706 2386526483 : ph3d(2,ig,iatblk)=ph12r*ph3i+ph12i*ph3r
707 : end do
708 : !$OMP END PARALLEL DO
709 : end do
710 :
711 8172161 : ABI_FREE(ph1kxred)
712 :
713 8172161 : end subroutine ph1d3d
714 : !!***
715 :
716 : !!****f* m_kg/getph
717 : !! NAME
718 : !! getph
719 : !!
720 : !! FUNCTION
721 : !! Compute three factors of one-dimensional structure factor phase
722 : !! for input atomic coordinates, for all planewaves which fit in fft box.
723 : !! The storage of these atomic factors is made according to the
724 : !! values provided by the index table atindx. This will save time in nonlop.
725 : !!
726 : !! INPUTS
727 : !! atindx(natom)=index table for atoms (see gstate.f)
728 : !! natom=number of atoms in cell.
729 : !! n1,n2,n3=dimensions of fft box (ngfft(3)).
730 : !! xred(3,natom)=reduced atomic coordinates.
731 : !!
732 : !! OUTPUT
733 : !! ph1d(2,(2*n1+1)*natom+(2*n2+1)*natom+(2*n3+1)*natom)=exp(2Pi i G.xred) for
734 : !! integer vector G with components ranging from -nj <= G <= nj.
735 : !! Real and imag given in usual Fortran convention.
736 : !!
737 : !! SOURCE
738 :
739 40153 : subroutine getph(atindx, natom, n1, n2, n3, ph1d, xred)
740 :
741 : !Arguments ------------------------------------
742 : !scalars
743 : integer,intent(in) :: n1,n2,n3,natom
744 : !arrays
745 : integer,intent(in) :: atindx(natom)
746 : real(dp),intent(in) :: xred(3,natom)
747 : real(dp),intent(out) :: ph1d(:,:)
748 :
749 : !Local variables-------------------------------
750 : !scalars
751 : integer,parameter :: im=2,re=1
752 : integer :: i1,i2,i3,ia,ii,ph1d_size1,ph1d_size2,ph1d_sizemin
753 : !character(len=500) :: msg
754 : real(dp) :: arg
755 : ! *************************************************************************
756 :
757 40153 : ph1d_size1=size(ph1d,1); ph1d_size2=size(ph1d,2)
758 40153 : ph1d_sizemin = (2*n1+1+2*n2+1+2*n3+1)*natom
759 40153 : if (ph1d_size1 /= 2 .or. ph1d_size2 < ph1d_sizemin) then
760 0 : ABI_BUG('Wrong ph1d sizes!')
761 : end if
762 :
763 117292 : do ia=1,natom
764 :
765 77139 : if(atindx(ia)<1 .or. natom<atindx(ia))then
766 0 : ABI_BUG('Wrong atindx(ia)!')
767 : endif
768 :
769 : ! Store the phase factor of atom number ia in place atindx(ia)
770 77139 : i1=(atindx(ia)-1)*(2*n1+1)
771 77139 : i2=(atindx(ia)-1)*(2*n2+1)+natom*(2*n1+1)
772 77139 : i3=(atindx(ia)-1)*(2*n3+1)+natom*(2*n1+1+2*n2+1)
773 :
774 3230196 : do ii=1,2*n1+1
775 3153057 : arg=two_pi*dble(ii-1-n1)*xred(1,ia)
776 3153057 : ph1d(re,ii+i1)=dcos(arg)
777 3230196 : ph1d(im,ii+i1)=dsin(arg)
778 : end do
779 :
780 3131674 : do ii=1,2*n2+1
781 3054535 : arg=two_pi*dble(ii-1-n2)*xred(2,ia)
782 3054535 : ph1d(re,ii+i2)=dcos(arg)
783 3131674 : ph1d(im,ii+i2)=dsin(arg)
784 : end do
785 :
786 3290185 : do ii=1,2*n3+1
787 3172893 : arg=two_pi*dble(ii-1-n3)*xred(3,ia)
788 3172893 : ph1d(re,ii+i3)=dcos(arg)
789 3250032 : ph1d(im,ii+i3)=dsin(arg)
790 : end do
791 :
792 : end do
793 :
794 : ! This to avoid uninitialized ph1d values
795 21727909 : if (ph1d_sizemin < ph1d_size2) ph1d(:,ph1d_sizemin+1:ph1d_size2)=zero
796 :
797 40153 : end subroutine getph
798 : !!***
799 :
800 : !!****f* m_kg/kpgstr
801 : !! NAME
802 : !! kpgstr
803 : !!
804 : !! FUNCTION
805 : !! Compute elements of the derivative the kinetic energy operator in reciprocal
806 : !! space at given k point wrt a single cartesian strain component
807 : !!
808 : !! INPUTS
809 : !! ecut=cut-off energy for plane wave basis sphere (Ha)
810 : !! ecutsm=smearing energy for plane wave kinetic energy (Ha)
811 : !! effmass_free=effective mass for electrons (1. in common case)
812 : !! gmet(3,3) = reciprocal lattice metric tensor (Bohr**-2)
813 : !! gprimd(3,3)=reciprocal space dimensional primitive translations
814 : !! istr=1,...6 specifies cartesian strain component 11,22,33,32,31,21
815 : !! kg(3,npw) = integer coordinates of planewaves in basis sphere.
816 : !! kpt(3) = reduced coordinates of k point
817 : !! npw = number of plane waves at kpt.
818 : !!
819 : !! OUTPUT
820 : !! dkinpw(npw)=d/deps(istr) ( (1/2)*(2 pi)**2 * (k+G)**2 )
821 : !!
822 : !! NOTES
823 : !! Src_6response/kpg3.f
824 : !!
825 : !! SOURCE
826 :
827 363703 : subroutine kpgstr(dkinpw,ecut,ecutsm,effmass_free,gmet,gprimd,istr,kg,kpt,npw)
828 :
829 : !Arguments -------------------------------
830 : !scalars
831 : integer,intent(in) :: istr,npw
832 : real(dp),intent(in) :: ecut,ecutsm,effmass_free
833 : !arrays
834 : integer,intent(in) :: kg(3,npw)
835 : real(dp),intent(in) :: gmet(3,3),gprimd(3,3),kpt(3)
836 : real(dp),intent(out) :: dkinpw(npw)
837 :
838 : !Local variables -------------------------
839 : !scalars
840 : integer :: ig,ii,ka,kb
841 : real(dp) :: dfsm,dkinetic,dkpg2,ecutsm_inv,fsm,gpk1,gpk2,gpk3,htpisq
842 : ! real(dp) :: d2fsm ! used in commented section below
843 : real(dp) :: kpg2,xx
844 : character(len=500) :: msg
845 : !arrays
846 : integer,save :: idx(12)=(/1,1,2,2,3,3,3,2,3,1,2,1/)
847 : real(dp) :: dgmetds(3,3)
848 :
849 : ! *********************************************************************
850 :
851 : !htpisq is (1/2) (2 Pi) **2:
852 363703 : htpisq=0.5_dp*(two_pi)**2
853 :
854 363703 : ecutsm_inv=0.0_dp
855 363703 : if(ecutsm>1.0d-20)ecutsm_inv=1/ecutsm
856 :
857 : !Compute derivative of metric tensor wrt strain component istr
858 363703 : if(istr<1 .or. istr>6)then
859 : write(msg, '(a,i10,a,a,a)' )&
860 0 : 'Input istr=',istr,' not allowed.',ch10,&
861 0 : 'Possible values are 1,2,3,4,5,6 only.'
862 0 : ABI_BUG(msg)
863 : end if
864 :
865 363703 : ka=idx(2*istr-1);kb=idx(2*istr)
866 1454812 : do ii = 1,3
867 4728139 : dgmetds(:,ii)=-(gprimd(ka,:)*gprimd(kb,ii)+gprimd(kb,:)*gprimd(ka,ii))
868 : end do
869 : !For historical reasons:
870 4728139 : dgmetds(:,:)=0.5_dp*dgmetds(:,:)
871 :
872 39982000 : do ig=1,npw
873 39618297 : gpk1=dble(kg(1,ig))+kpt(1)
874 39618297 : gpk2=dble(kg(2,ig))+kpt(2)
875 39618297 : gpk3=dble(kg(3,ig))+kpt(3)
876 : kpg2=htpisq*&
877 : & ( gmet(1,1)*gpk1**2+ &
878 : & gmet(2,2)*gpk2**2+ &
879 : & gmet(3,3)*gpk3**2 &
880 : & +2.0_dp*(gpk1*gmet(1,2)*gpk2+ &
881 : & gpk1*gmet(1,3)*gpk3+ &
882 39618297 : & gpk2*gmet(2,3)*gpk3 ) )
883 : dkpg2=htpisq*2.0_dp*&
884 : & (gpk1*(dgmetds(1,1)*gpk1+dgmetds(1,2)*gpk2+dgmetds(1,3)*gpk3)+ &
885 : & gpk2*(dgmetds(2,1)*gpk1+dgmetds(2,2)*gpk2+dgmetds(2,3)*gpk3)+ &
886 39618297 : & gpk3*(dgmetds(3,1)*gpk1+dgmetds(3,2)*gpk2+dgmetds(3,3)*gpk3) )
887 39618297 : dkinetic=dkpg2
888 39618297 : if(kpg2>ecut-ecutsm)then
889 4693893 : if(kpg2>ecut-tol12)then
890 : ! The wavefunction has been filtered : no derivative
891 : dkinetic=0.0_dp
892 : else
893 4693893 : xx=(ecut-kpg2)*ecutsm_inv
894 : ! This kinetic cutoff smoothing function and its xx derivatives
895 : ! were produced with Mathematica and the fortran code has been
896 : ! numerically checked against Mathematica.
897 4693893 : fsm=1.0_dp/(xx**2*(3+xx*(1+xx*(-6+3*xx))))
898 4693893 : dfsm=-3.0_dp*(-1+xx)**2*xx*(2+5*xx)*fsm**2
899 : ! d2fsm=6.0_dp*xx**2*(9+xx*(8+xx*(-52+xx*(-3+xx*(137+xx*&
900 : ! & (-144+45*xx))))))*fsm**3
901 4693893 : dkinetic=dkpg2*(fsm-ecutsm_inv*kpg2*dfsm)
902 : end if
903 : end if
904 39982000 : dkinpw(ig)=dkinetic/effmass_free
905 : end do
906 :
907 363703 : end subroutine kpgstr
908 : !!***
909 :
910 : !!****f* m_kg/mkkpg
911 : !! NAME
912 : !! mkkpg
913 : !!
914 : !! FUNCTION
915 : !! Compute all (k+G) vectors (dp, in reduced coordinates) for given k point,
916 : !! from integer coordinates of G vectors.
917 : !! Eventually compute related data.
918 : !!
919 : !! INPUTS
920 : !! kg(3,npw)=integer coords of planewaves in basis sphere
921 : !! kpt(3)=k point in terms of recip. translations
922 : !! nkpg=second dimension of array kpg
923 : !! npw=number of plane waves in reciprocal space
924 : !!
925 : !! OUTPUT
926 : !! kpg(npw,3)= (k+G) components
927 : !! === if nkpg==9 ===
928 : !! kpg(npw,4:9)= [(k+G)_a].[(k+G)_b] quantities
929 : !!
930 : !! SOURCE
931 :
932 :
933 3419465 : subroutine mkkpg(kg, kpg, kpt, nkpg, npw)
934 :
935 : !Arguments ------------------------------------
936 : !scalars
937 : integer,intent(in) :: nkpg,npw
938 : !arrays
939 : integer,intent(in) :: kg(3,npw)
940 : real(dp),intent(in) :: kpt(3)
941 : real(dp),intent(out) :: kpg(npw,nkpg)
942 :
943 : !Local variables-------------------------------
944 : !scalars
945 : integer :: ipw,mu,mua,mub
946 : character(len=500) :: msg
947 : !arrays
948 : integer,parameter :: alpha(6)=(/1,2,3,3,3,2/),beta(6)=(/1,2,3,2,1,1/)
949 : ! *************************************************************************
950 :
951 3419465 : if (nkpg==0) return
952 :
953 : !-- Test nkpg --
954 3419369 : if (nkpg/=3.and.nkpg/=9) then
955 0 : write(msg, '(a,i0)' )' Bad value for nkpg !',nkpg
956 0 : ABI_BUG(msg)
957 : end if
958 :
959 : !-- Compute (k+G) --
960 : !$OMP PARALLEL DO COLLAPSE(2) &
961 : !$OMP PRIVATE(mu,ipw)
962 418056868 : do ipw=1,npw
963 1661969365 : do mu=1,3
964 1658549996 : kpg(ipw,mu)=kpt(mu)+dble(kg(mu,ipw))
965 : end do
966 : end do
967 : !$OMP END PARALLEL DO
968 :
969 : !-- Compute [(k+G)_a].[(k+G)_b] --
970 3419369 : if (nkpg==9) then
971 : !$OMP PARALLEL DO COLLAPSE(2) &
972 : !$OMP PRIVATE(ipw,mu,mua,mub)
973 2884396 : do ipw=1,npw
974 20047636 : do mu=4,9
975 17163240 : mua=alpha(mu-3);mub=beta(mu-3)
976 20023780 : kpg(ipw,mu)=kpg(ipw,mua)*kpg(ipw,mub)
977 : end do
978 : end do
979 : !$OMP END PARALLEL DO
980 : end if
981 :
982 : end subroutine mkkpg
983 : !!***
984 :
985 : !!****f* ABINIT/mkpwind_k
986 : !! NAME
987 : !! mkpwind_k
988 : !!
989 : !! FUNCTION
990 : !! Make plane wave index at k point for basis at second k point,
991 : !! needed to compute overlaps $\langle u_{k,n}|u_{k+b,n}\rangle$
992 : !! as appear in Berry phase derived quantities
993 : !!
994 : !! INPUTS
995 : !! dk(3)=real vector difference of ket kpt - bra kpt
996 : !! dtset <type(dataset_type)>=all input variables in this dataset
997 : !! fnkpt=number of kpts in full BZ
998 : !! fkptns=kpts in full BZ
999 : !! gmet(3,3)=metric in reciprocal space
1000 : !! indkk_f2ibz(fnkpt,6)=information on folding from FBZ to IBZ (see initberry or initorbmag)
1001 : !! ikpt=index of bra k pt in FBZ
1002 : !! ikpt1=index of neighbour ket k pt in FBZ
1003 : !! mpi_enreg=information about MPI parallelization
1004 : !! npwarr(dtset%nkpt)=npw at each kpt
1005 : !! symrec(3,3,nsym) = symmetries in reciprocal space in terms of
1006 : !! reciprocal space primitive translations
1007 : !!
1008 : !! OUTPUT
1009 : !! pwind_k1(dtset%mpw)=output index of ikpt1 basis states refered to ikpt
1010 : !!
1011 : !! SOURCE
1012 :
1013 0 : subroutine mkpwind_k(dk,dtset,fnkpt,fkptns,gmet,indkk_f2ibz,ikpt,ikpt1,&
1014 0 : & mpi_enreg,npwarr,pwind_k1,symrec)
1015 :
1016 : !Arguments ------------------------------------
1017 : !scalars
1018 : integer,intent(in) :: fnkpt,ikpt,ikpt1
1019 : type(dataset_type),intent(in) :: dtset
1020 : type(MPI_type), intent(inout) :: mpi_enreg
1021 :
1022 : !arrays
1023 : integer,intent(in) :: indkk_f2ibz(fnkpt,6)
1024 : integer,intent(in) :: npwarr(dtset%nkpt)
1025 : integer,intent(in) :: symrec(3,3,dtset%nsym)
1026 : integer,intent(out) :: pwind_k1(dtset%mpw)
1027 : real(dp),intent(in) :: dk(3),fkptns(3,fnkpt),gmet(3,3)
1028 :
1029 : !Local variables -------------------------
1030 : !scalars
1031 : integer :: exchn2n3d,idum1,ikg1,ikpti,ikpt1i,ipw,istwf_k,isym,isym1,jpw,npw_k,npw_k1
1032 : real(dp) :: ecut_eff
1033 :
1034 : !arrays
1035 0 : integer,allocatable :: kg_k(:,:),kg1_k(:,:)
1036 : real(dp) :: dg(3),dum33(3,3),kpt(3),kpt1(3),iadum(3),iadum1(3)
1037 : ! ***********************************************************************
1038 :
1039 0 : ikpti = indkk_f2ibz(ikpt,1)
1040 0 : ikpt1i = indkk_f2ibz(ikpt1,1)
1041 :
1042 0 : ecut_eff = dtset%ecut*(dtset%dilatmx)**2
1043 0 : exchn2n3d = 0 ; istwf_k = 1 ; ikg1 = 0
1044 :
1045 : ! Build basis sphere of plane waves for the k-point
1046 : ! we avoid using the global kg data because of difficulties in parallel-ism
1047 0 : ABI_MALLOC(kg_k,(3,dtset%mpw))
1048 0 : kg_k(:,:) = 0
1049 0 : kpt(:) = dtset%kptns(:,ikpti)
1050 0 : call kpgsph(ecut_eff,exchn2n3d,gmet,ikg1,ikpt,istwf_k,kg_k,kpt,1,mpi_enreg,dtset%mpw,npw_k)
1051 :
1052 : ! Build basis sphere of plane waves for the nearest neighbour of the k-point
1053 0 : ABI_MALLOC(kg1_k,(3,dtset%mpw))
1054 0 : kg1_k(:,:) = 0
1055 0 : kpt1(:) = dtset%kptns(:,ikpt1i)
1056 0 : call kpgsph(ecut_eff,exchn2n3d,gmet,ikg1,ikpt,istwf_k,kg1_k,kpt1,1,mpi_enreg,dtset%mpw,npw_k1)
1057 :
1058 : ! Deal with symmetry transformations
1059 :
1060 : ! bra k-point k(b) and IBZ k-point kIBZ(b) related by
1061 : ! k(b) = alpha(b) S(b)^t kIBZ(b) + G(b)
1062 : ! where alpha(b), S(b) and G(b) are given by indkk_f2ibz
1063 : !
1064 : ! For the ket k-point:
1065 : ! k(k) = alpha(k) S(k)^t kIBZ(k) + G(k) - GBZ(k)
1066 : ! where GBZ(k) takes k(k) to the BZ
1067 :
1068 0 : isym = indkk_f2ibz(ikpt,2)
1069 0 : isym1 = indkk_f2ibz(ikpt1,2)
1070 :
1071 : ! Construct transformed G vector that enters the matching condition:
1072 : ! alpha(k) S(k)^{t,-1} ( -G(b) - GBZ(k) + G(k) )
1073 :
1074 : dg(:) = -indkk_f2ibz(ikpt,3:5) &
1075 : & - nint(-fkptns(:,ikpt) - dk(:) - tol10 + fkptns(:,ikpt1)) &
1076 0 : & + indkk_f2ibz(ikpt1,3:5)
1077 :
1078 0 : iadum(:) = MATMUL(TRANSPOSE(dtset%symrel(:,:,isym1)),dg(:))
1079 :
1080 0 : dg(:) = iadum(:)
1081 :
1082 : ! Construct S(k)^{t,-1} S(b)^{t}
1083 :
1084 0 : dum33(:,:) = MATMUL(TRANSPOSE(dtset%symrel(:,:,isym1)),symrec(:,:,isym))
1085 :
1086 : ! Construct alpha(k) alpha(b)
1087 :
1088 0 : pwind_k1(:) = 0
1089 0 : npw_k = npwarr(ikpti)
1090 0 : do ipw = 1, npw_k
1091 :
1092 : ! NOTE: the bra G vector is taken for the sym-related IBZ k point,
1093 : ! not for the FBZ k point
1094 :
1095 : ! original code from initberry
1096 : ! iadum(:) = kg(:,kgindex(ikpti) + ipw)
1097 :
1098 0 : iadum(:) = kg_k(:,ipw)
1099 :
1100 : ! to determine r.l.v. matchings, we transformed the bra vector
1101 : ! Rotation
1102 0 : iadum1(:)=0
1103 0 : do idum1=1,3
1104 0 : iadum1(:)=iadum1(:)+dum33(:,idum1)*iadum(idum1)
1105 : end do
1106 0 : iadum(:)=iadum1(:)
1107 0 : iadum(:) = iadum(:) + dg(:)
1108 :
1109 0 : do jpw = 1, npw_k1
1110 0 : iadum1(1:3) = kg1_k(1:3,jpw)
1111 : if ( (iadum(1) == iadum1(1)).and. &
1112 0 : & (iadum(2) == iadum1(2)).and. &
1113 0 : & (iadum(3) == iadum1(3)) ) then
1114 0 : pwind_k1(ipw) = jpw
1115 0 : exit
1116 : end if
1117 : end do
1118 : end do
1119 :
1120 0 : ABI_FREE(kg_k)
1121 0 : ABI_FREE(kg1_k)
1122 :
1123 0 : end subroutine mkpwind_k
1124 : !!***
1125 :
1126 : !!****f* m_kg/mkkpgcart
1127 : !! NAME
1128 : !! mkkpgcart
1129 : !!
1130 : !! FUNCTION
1131 : !! Compute all (k+G) vectors (dp, in cartesian coordinates) for given k point,
1132 : !! from integer coordinates of G vectors.
1133 : !!
1134 : !! INPUTS
1135 : !! gprimd(3,3)=dimensional reciprocal space primitive translations
1136 : !! kg(3,npw)=integer coords of planewaves in basis sphere
1137 : !! kpt(3)=k point in terms of recip. translations
1138 : !! nkpg=second dimension of array kpg
1139 : !! npw=number of plane waves in reciprocal space
1140 : !!
1141 : !! OUTPUT
1142 : !! kpg(npw,3)= (k+G) components
1143 : !!
1144 : !! SOURCE
1145 :
1146 746496 : subroutine mkkpgcart(gprimd,kg,kpgcar,kpt,nkpg,npw)
1147 :
1148 : !Arguments ------------------------------------
1149 : !scalars
1150 : integer,intent(in) :: nkpg,npw
1151 : !arrays
1152 : real(dp),intent(in) :: gprimd(3,3)
1153 : integer,intent(in) :: kg(3,npw)
1154 : real(dp),intent(in) :: kpt(3)
1155 : real(dp),intent(out) :: kpgcar(npw,nkpg)
1156 :
1157 : !Local variables-------------------------------
1158 : !scalars
1159 : integer :: ipw,mu
1160 : character(len=500) :: msg
1161 : !arrays
1162 746496 : real(dp),allocatable :: kpg(:,:)
1163 : ! *************************************************************************
1164 :
1165 : DBG_ENTER("COLL")
1166 :
1167 746496 : if (nkpg==0) return
1168 :
1169 : !-- Test nkpg --
1170 746496 : if (nkpg/=3) then
1171 0 : write(msg, '(a,i0)' )' Bad value for nkpg !',nkpg
1172 0 : ABI_BUG(msg)
1173 : end if
1174 :
1175 : !-- Compute (k+G) --
1176 2985984 : ABI_MALLOC(kpg,(npw,nkpg))
1177 : !$OMP PARALLEL DO COLLAPSE(2)
1178 74486304 : do ipw=1,npw
1179 295705728 : do mu=1,3
1180 294959232 : kpg(ipw,mu)=kpt(mu)+dble(kg(mu,ipw))
1181 : end do
1182 : end do
1183 : !$OMP END PARALLEL DO
1184 :
1185 : !$OMP PARALLEL DO
1186 74486304 : do ipw=1,npw
1187 73739808 : kpgcar(ipw,1)=kpg(ipw,1)*gprimd(1,1)+kpg(ipw,2)*gprimd(1,2)+kpg(ipw,3)*gprimd(1,3)
1188 73739808 : kpgcar(ipw,2)=kpg(ipw,1)*gprimd(2,1)+kpg(ipw,2)*gprimd(2,2)+kpg(ipw,3)*gprimd(2,3)
1189 74486304 : kpgcar(ipw,3)=kpg(ipw,1)*gprimd(3,1)+kpg(ipw,2)*gprimd(3,2)+kpg(ipw,3)*gprimd(3,3)
1190 : end do
1191 : !$OMP END PARALLEL DO
1192 :
1193 746496 : ABI_FREE(kpg)
1194 : DBG_EXIT("COLL")
1195 :
1196 : end subroutine mkkpgcart
1197 : !!***
1198 :
1199 : !!****f* ABINIT/mkkin_metdqdq
1200 : !! NAME
1201 : !! mkkin_metdqdq
1202 : !!
1203 : !! FUNCTION
1204 : !! Compute elements of the second q-gradient of the metric
1205 : !! kinetic energy operator in reciprocal
1206 : !! space at given k point wrt cartesian q components.
1207 : !!
1208 : !! INPUTS
1209 : !! effmass=effective mass for electrons (1. in common case)
1210 : !! gprimd(3,3)=reciprocal space dimensional primitive translations
1211 : !! idir= strain perturbation direction
1212 : !! kg(3,npw) = integer coordinates of planewaves in basis sphere.
1213 : !! kpt(3) = reduced coordinates of k point
1214 : !! npw = number of plane waves at kpt.
1215 : !! qdir = direction of the first q-gradient
1216 : !!
1217 : !! OUTPUT
1218 : !! dqdqkinpw(npw)=d/deps(istr) ( (1/2)*(2 pi)**2 * (k+G)**2 )
1219 : !!
1220 : !! NOTES
1221 : !! **Since the 2nd derivative w.r.t q-vector is calculated along cartesian
1222 : !! directions, the 1/twopi**2 factor (that in the rest of the code is applied
1223 : !! in the reduced to cartesian derivative conversion process) is here
1224 : !! explicictly included in the formulas.
1225 : !!
1226 : !! **Notice that idir=1-9, in contrast to the strain perturbation (idir=1-6),
1227 : !! because this term is not symmetric w.r.t permutations of the two strain
1228 : !! indices.
1229 : !!
1230 : !! **A -i factor has been factorized out in all the contributions of the second
1231 : !! q-gradient of the metric Hamiltonian. This is lately included in the contribution
1232 : !! of the corresponing term (T4) to the flexoelectric tensor in dfpt_flexoout.F90
1233 : !!
1234 : !! SOURCE
1235 :
1236 93312 : subroutine mkkin_metdqdq(dqdqkinpw,effmass,gprimd,idir,kg,kpt,npw,qdir)
1237 :
1238 : !Arguments -------------------------------
1239 : !scalars
1240 : integer,intent(in) :: idir,npw,qdir
1241 : real(dp),intent(in) :: effmass
1242 : !arrays
1243 : integer,intent(in) :: kg(3,npw)
1244 : real(dp),intent(in) :: gprimd(3,3),kpt(3)
1245 : real(dp),intent(out) :: dqdqkinpw(npw)
1246 :
1247 : !Local variables -------------------------
1248 : !scalars
1249 : integer :: beta,delta,gamma,ig,ka,kb
1250 : real(dp) :: delbd,delbg,deldg, dkinetic,gpk1,gpk2,gpk3,htpi
1251 : !arrays
1252 : real(dp) :: gpkc(3)
1253 : integer,parameter :: idx(18) = [1,1,2,2,3,3,3,2,3,1,2,1,2,3,1,3,1,2]
1254 : ! *********************************************************************
1255 :
1256 : !htpi is (1/2) (2 Pi):
1257 93312 : htpi=0.5_dp*two_pi
1258 :
1259 93312 : ka=idx(2*idir-1);kb=idx(2*idir)
1260 :
1261 : !For easier formula implementation
1262 93312 : beta=ka
1263 93312 : delta=kb
1264 93312 : gamma=qdir
1265 :
1266 : !Kronecker deltas
1267 93312 : delbd=0.0_dp; delbg=0.0_dp; deldg=0.0_dp
1268 93312 : if (beta==delta) delbd=1.0_dp
1269 93312 : if (beta==gamma) delbg=1.0_dp
1270 93312 : if (delta==gamma) deldg=1.0_dp
1271 :
1272 9310788 : do ig=1,npw
1273 9217476 : gpk1=dble(kg(1,ig))+kpt(1)
1274 9217476 : gpk2=dble(kg(2,ig))+kpt(2)
1275 9217476 : gpk3=dble(kg(3,ig))+kpt(3)
1276 :
1277 : ! Obtain G in cartesian coordinates
1278 9217476 : gpkc(1)=gprimd(1,1)*gpk1+gprimd(1,2)*gpk2+gprimd(1,3)*gpk3
1279 9217476 : gpkc(2)=gprimd(2,1)*gpk1+gprimd(2,2)*gpk2+gprimd(2,3)*gpk3
1280 9217476 : gpkc(3)=gprimd(3,1)*gpk1+gprimd(3,2)*gpk2+gprimd(3,3)*gpk3
1281 :
1282 9217476 : dkinetic=htpi*(2.0_dp*deldg*gpkc(beta)+delbg*gpkc(delta)+ delbd*gpkc(gamma))
1283 :
1284 9310788 : dqdqkinpw(ig)=dkinetic/effmass
1285 : end do
1286 :
1287 93312 : end subroutine mkkin_metdqdq
1288 : !!***
1289 :
1290 0 : end module m_kg
1291 : !!***
|
