Branch data 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 : 38746 : 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 : 19373 : ecut_pw=ecut
119 [ - + ]: 19373 : 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 : 19373 : call bound(largesq,boxsq,gbound,gmet,kpt,ngfft,plane)
126 : 19373 : effcut=0.5_dp * (two_pi)**2 * boxsq
127 : 19373 : sphsq=2._dp*ecut_pw/two_pi**2
128 : :
129 [ + + ]: 19373 : 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 : 19370 : boxcut = sqrt(boxsq/sphsq)
146 : 19370 : cutrad = min(2.0_dp,boxcut)
147 : 19370 : gsqcut = (cutrad**2)*(2.0_dp*ecut_pw)/two_pi**2
148 : :
149 [ + - ]: 19370 : if(ecut>-tol8)then
150 : :
151 : 19370 : write(msg, '(a,a,3f8.4,a,3i4,a,a,f11.3,3x,a,f10.5)' ) ch10,&
152 : 19370 : ' getcut: wavevector=',kpt,' ngfft=',ngfft(1:3),ch10,&
153 : 38740 : ' ecut(hartree)=',ecut+tol8,'=> boxcut(ratio)=',boxcut+tol8
154 [ + + ]: 19370 : if(iout/=std_out) call wrtout(iout,msg)
155 : 19370 : call wrtout(std_out,msg)
156 : :
157 [ - + ]: 19370 : 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 [ + + ]: 19370 : if (boxcut>2.2_dp) then
172 : 7196 : write(msg, '(a,a,a,a,a,a,a,a,a,a,a,f12.6,a,a)' ) ch10,&
173 : 7196 : ' getcut : COMMENT -',ch10,&
174 : 7196 : ' Note that boxcut > 2.2 ; recall that',' boxcut=Gcut(box)/Gcut(sphere) = 2',ch10,&
175 : 7196 : ' is sufficient for exact treatment of convolution.',ch10,&
176 : 7196 : ' Such a large boxcut is a waste : you could raise ecut',ch10,&
177 : 14392 : ' e.g. ecut=',effcut*0.25_dp+tol8,' Hartrees makes boxcut=2',ch10
178 [ + + ]: 7196 : if(iout/=std_out) call wrtout(iout,msg)
179 : 7196 : call wrtout(std_out,msg)
180 : : end if
181 : :
182 [ + + ]: 19370 : 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 : 19373 : 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 : 12578 : 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 : 12578 : 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 [ - + ]: 12578 : ABI_MALLOC(kg,(3,100))
259 : :
260 : : !set mpw to zero, as needed for only counting in kpgsph
261 : 12578 : mpw = 0
262 : :
263 : : !Might be parallelized over k points ? !
264 [ + + ]: 477325 : do ikpt = 1,nkpt
265 : : ! Do computation of G sphere, returning npw
266 [ + + ]: 1858988 : kpoint(:)=kptns(:,ikpt)
267 : 464747 : istwf_k=istwfk(ikpt)
268 : 464747 : call kpgsph(ecut,exchn2n3d,gmet,0,ikpt,istwf_k,kg,kpoint,0,mpi_enreg,0,npw)
269 : 477325 : mpw = max(npw,mpw)
270 : : end do
271 : :
272 : 12578 : write(msg,'(a,i0)') ' getmpw: optimal value of mpw= ',mpw
273 : 12578 : call wrtout(std_out,msg)
274 : :
275 [ - + ]: 12578 : ABI_FREE(kg)
276 : :
277 : 12578 : 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 : 3654240 : 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 : 3654240 : l_vecpot = .false.
350 [ + + ]: 3654240 : 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 : 3654240 : htpisq=0.5_dp*(two_pi)**2
356 : :
357 : 3654240 : ecutsm_inv=0.0_dp
358 [ + + ]: 3654240 : if(ecutsm>1.0d-20)ecutsm_inv=1/ecutsm
359 : :
360 : 3654240 : gmet_break(:,:)=gmet(:,:)
361 : 3654240 : gmet_break(1,1)=(1.0_dp+break_symm)*gmet(1,1)
362 : 3654240 : gmet_break(3,3)=(1.0_dp-break_symm)*gmet(3,3)
363 : :
364 : 3654240 : order=0 ! Compute the kinetic operator
365 [ + + ]: 3654240 : if (idir1>0.and.idir1<4) then
366 : 350750 : order=1 ! Compute the 1st derivative of the kinetic operator
367 [ + + ]: 350750 : 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 [ + + ]: 3654240 : 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 [ + + ]: 799926082 : do ig=1,npw
387 : 796271842 : gpk1=dble(kg(1,ig))+kpt(1)
388 : 796271842 : gpk2=dble(kg(2,ig))+kpt(2)
389 : 796271842 : 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 : 796271842 : & gpk2*gmet_break(2,3)*gpk3 ) )
397 [ + + ]: 796271842 : 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 : 746391007 : select case (order)
407 : : case (0)
408 : 746391007 : kinetic=kpg2
409 [ + + ]: 746391007 : if (l_vecpot) kinetic=kinetic+two_pi*akpg+0.5_dp*asq
410 : : case (1)
411 : : dkpg2=htpisq*2.0_dp*&
412 : 48263855 : & (gmet_break(idir1,1)*gpk1+gmet_break(idir1,2)*gpk2+gmet_break(idir1,3)*gpk3)
413 : 48263855 : kinetic=dkpg2
414 [ - + ]: 48263855 : 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 [ + + + ]: 796271842 : kinetic=dkpg2
420 : : end select
421 : :
422 [ + + ]: 796271842 : 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 : 799926082 : kinpw(ig)=kinetic/effmass_free
463 : : end do
464 : : !$OMP END PARALLEL DO
465 : :
466 : 3654240 : 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 : 11642 : subroutine kpgio(ecut,exchn2n3d,gmet,istwfk,kg,kptns,mkmem,nband,nkpt,&
505 : 11642 : & 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 : 11642 : me=mpi_enreg%me_kpt
529 : :
530 [ + + - + ]: 11642 : if((mpi_enreg%paralbd==1) .and. (mode_paral=='PERS')) then
531 [ + + ]: 5538 : 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 [ + + ]: 441487 : npwarr(:)=0
550 [ + + ]: 441487 : npwtot(:)=0
551 : :
552 [ + + + + ]: 225734830 : kg=0
553 : 11642 : ikg=0
554 : : !Find (k+G) sphere for each k.
555 : :
556 [ + + ]: 441487 : do ikpt=1,nkpt
557 : :
558 : 429845 : nband_k = nband(ikpt)
559 : :
560 [ + + ]: 429845 : if(mode_paral=='PERS')then
561 [ + + ]: 429767 : if(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,-1,me)) cycle
562 : : end if
563 : :
564 [ + + ]: 1406036 : kpoint(:)=kptns(:,ikpt)
565 : 351509 : istwf_k=istwfk(ikpt)
566 : 351509 : call kpgsph(ecut,exchn2n3d,gmet,ikg,ikpt,istwf_k,kg,kpoint,mkmem,mpi_enreg,mpw,npw1)
567 : :
568 : 351509 : test_npw=.true.
569 : : if (xmpi_paral==1)then
570 [ + + ]: 351509 : if (mode_paral=='PERS')then
571 [ + + + + ]: 2829261 : test_npw=(minval(mpi_enreg%proc_distrb(ikpt,1:nband_k,1:nsppol))==me)
572 : : end if
573 : : end if
574 [ + + ]: 351509 : 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 : 792996 : 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 [ + + ]: 11642 : if(mode_paral == 'PERS') then
594 : 11564 : 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 [ + + ]: 441487 : npwtot(:)=npwarr(:)
602 : :
603 : : !Taking into account istwfk
604 [ + + ]: 441487 : do ikpt=1,nkpt
605 [ + + ]: 441487 : if(istwfk(ikpt)>1)then
606 [ + + ]: 3155 : if(istwfk(ikpt)==2)then
607 : 1459 : 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 : 11642 : 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 : 8174365 : 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 : 8174365 : real(dp),allocatable :: ph1kxred(:,:)
662 : : ! *************************************************************************
663 : :
664 [ - + ]: 8174365 : 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 [ - + + - : 24523095 : ABI_MALLOC(ph1kxred,(2,-n1:n1))
+ - + - -
+ + - ]
674 : :
675 : : ! ia runs from iatom to jatom
676 [ + + ]: 18457335 : do ia=iatom,jatom
677 : :
678 : : ! iatblk runs from 1 to matblk
679 : 10282970 : iatblk=ia-iatom+1
680 : 10282970 : shift1=1+n1+(ia-1)*(2*n1+1)
681 : 10282970 : shift2=1+n2+(ia-1)*(2*n2+1)+natom*(2*n1+1)
682 : 10282970 : 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 : 10282970 : phkxr=phkxred(1,ia)
685 : 10282970 : phkxi=phkxred(2,ia)
686 [ + + ]: 338598222 : do i1=-n1,n1
687 : 328315252 : ph1kxred(1,i1)=ph1d(1,i1+shift1)*phkxr-ph1d(2,i1+shift1)*phkxi
688 : 338598222 : 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 [ + + ]: 2397503666 : do ig=1,npw_k
694 : 2379046331 : kg1=kg_k(1,ig)
695 : 2379046331 : kg2=kg_k(2,ig)+shift2
696 : 2379046331 : kg3=kg_k(3,ig)+shift3
697 : 2379046331 : ph1r=ph1kxred(1,kg1)
698 : 2379046331 : ph1i=ph1kxred(2,kg1)
699 : 2379046331 : ph2r=ph1d(1,kg2)
700 : 2379046331 : ph2i=ph1d(2,kg2)
701 : 2379046331 : ph3r=ph1d(1,kg3)
702 : 2379046331 : ph3i=ph1d(2,kg3)
703 : 2379046331 : ph12r=ph1r*ph2r-ph1i*ph2i
704 : 2379046331 : ph12i=ph1r*ph2i+ph1i*ph2r
705 : 2379046331 : ph3d(1,ig,iatblk)=ph12r*ph3r-ph12i*ph3i
706 : 2389329301 : ph3d(2,ig,iatblk)=ph12r*ph3i+ph12i*ph3r
707 : : end do
708 : : !$OMP END PARALLEL DO
709 : : end do
710 : :
711 [ - + ]: 8174365 : ABI_FREE(ph1kxred)
712 : :
713 : 8174365 : 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 [ - + ]: 40240 : 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 : 40240 : ph1d_size1=size(ph1d,1); ph1d_size2=size(ph1d,2)
758 : 40240 : ph1d_sizemin = (2*n1+1+2*n2+1+2*n3+1)*natom
759 [ - + ]: 40240 : if (ph1d_size1 /= 2 .or. ph1d_size2 < ph1d_sizemin) then
760 : 0 : ABI_BUG('Wrong ph1d sizes!')
761 : : end if
762 : :
763 [ + + ]: 117517 : do ia=1,natom
764 : :
765 [ + - - + ]: 77277 : 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 : 77277 : i1=(atindx(ia)-1)*(2*n1+1)
771 : 77277 : i2=(atindx(ia)-1)*(2*n2+1)+natom*(2*n1+1)
772 : 77277 : i3=(atindx(ia)-1)*(2*n3+1)+natom*(2*n1+1+2*n2+1)
773 : :
774 [ + + ]: 3240904 : do ii=1,2*n1+1
775 : 3163627 : arg=two_pi*dble(ii-1-n1)*xred(1,ia)
776 : 3163627 : ph1d(re,ii+i1)=dcos(arg)
777 : 3240904 : ph1d(im,ii+i1)=dsin(arg)
778 : : end do
779 : :
780 [ + + ]: 3142382 : do ii=1,2*n2+1
781 : 3065105 : arg=two_pi*dble(ii-1-n2)*xred(2,ia)
782 : 3065105 : ph1d(re,ii+i2)=dcos(arg)
783 : 3142382 : ph1d(im,ii+i2)=dsin(arg)
784 : : end do
785 : :
786 [ + + ]: 3300980 : do ii=1,2*n3+1
787 : 3183463 : arg=two_pi*dble(ii-1-n3)*xred(3,ia)
788 : 3183463 : ph1d(re,ii+i3)=dcos(arg)
789 : 3260740 : ph1d(im,ii+i3)=dsin(arg)
790 : : end do
791 : :
792 : : end do
793 : :
794 : : ! This to avoid uninitialized ph1d values
795 [ + + + + : 21727996 : if (ph1d_sizemin < ph1d_size2) ph1d(:,ph1d_sizemin+1:ph1d_size2)=zero
+ + ]
796 : :
797 : 40240 : 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 : 3421577 : 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 [ + + ]: 3421577 : if (nkpg==0) return
952 : :
953 : : !-- Test nkpg --
954 [ - + ]: 3421481 : 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 [ + + ]: 418684074 : do ipw=1,npw
963 [ + + ]: 1664471853 : do mu=1,3
964 : 1661050372 : 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 [ + + ]: 3421481 : 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 : : !!***
|
