Line data Source code
1 : !!****m* ABINIT/m_cgtk
2 : !! NAME
3 : !! m_cgtk
4 : !!
5 : !! FUNCTION
6 : !!
7 : !!
8 : !! COPYRIGHT
9 : !! Copyright (C) 2008-2026 ABINIT group (MG)
10 : !! This file is distributed under the terms of the
11 : !! GNU General Public License, see ~abinit/COPYING
12 : !! or http://www.gnu.org/copyleft/gpl.txt .
13 : !!
14 : !! SOURCE
15 :
16 : #if defined HAVE_CONFIG_H
17 : #include "config.h"
18 : #endif
19 :
20 : #include "abi_common.h"
21 :
22 : module m_cgtk
23 :
24 : use defs_basis
25 : use m_errors
26 : use m_abicore
27 : use m_xmpi
28 : use m_time
29 :
30 : use m_fstrings, only : itoa, sjoin
31 : use defs_abitypes, only : MPI_type
32 : use m_matrix, only : mati3inv
33 : use m_geometry, only : getspinrot
34 : use m_crystal, only : crystal_t
35 : use m_fftcore, only : sphere
36 : use m_kg, only : ph1d3d, getph
37 : use m_pawcprj, only : pawcprj_type, pawcprj_zaxpby
38 :
39 : implicit none
40 :
41 : private
42 : !!***
43 :
44 : !interface cgtk_rotate
45 : ! module procedure cgtk_rotate_dp
46 : ! module procedure cgtk_rotate_sp
47 : !end interface cgtk_rotate
48 :
49 : public :: cgtk_rotate ! Recostruct wfs in the BZ from IBZ using symmetry tables generated
50 : ! with the the symrel^T conventions.
51 : public :: cgtk_rotate_symrec ! Similar to cgtk_rotate but assumes symrec conventions.
52 : public :: cgtk_change_gsphere
53 : public :: cgtk_fixphase
54 : !!***
55 :
56 : integer,private,parameter :: to_box = 1, to_sph = -1, me_g0 = 1, ndat1 = 1
57 : integer,private,parameter :: no_shift(3) = 0
58 :
59 : contains
60 : !!***
61 :
62 : !!****f* ABINIT/cgtk_rotate
63 : !! NAME
64 : !! cgtk_rotate
65 : !!
66 : !! FUNCTION
67 : !! Reconstruct wavefunction cg2 in the BZ from the symmetrical image cg1 by applying a symmetry operation.
68 : !! Note that there are two possible conventions for mapping k-points:
69 : !!
70 : !! 1) k2 = T symrel(:,:, isym)^t k1 + g0 (note transpose of symrel)
71 : !!
72 : !! 2) k2 = T symrec(:,:, isym) k1 + g0
73 : !!
74 : !! where T is for time-reversal (itimrev)
75 : !!
76 : !! This routine assumes the FIRST convention that, unfortunately, is not very handy.
77 : !! The second convention, indeed, is the most natural one when mapping k-points.
78 : !!
79 : !! INPUTS
80 : !! cryst=crystalline structure
81 : !! kpt1(3)=k-point in cg1.
82 : !! isym=Index of symmetry operation (symrel^T convention)
83 : !! itimrev=1 if time-reversal is needed else 0.
84 : !! g0(3)=g0 vector
85 : !! nspinor=Number of spinor components.
86 : !! ndat=Number of wavefunctions
87 : !! npw1, npw2=Number of G-vectors in kg1 and kg2.
88 : !! kg1(3,npw1), kg2(3,npw2) = G vectors in cg1, and cg2.
89 : !! istwf1, istwf2= Storage mode for cg1 and cg2
90 : !! work_ngfft(18)= Specifies the size of the workspace array work.
91 : !! IMPORTANT: must be large enough to account for all possible shifts of the g-sphere.
92 : !! The caller is responsible for computing the max size needed to handle all the possible symmetrization.
93 : !! cg1(2, npw1, nspinor, ndat)=Wavefunctions in the IBZ
94 : !!
95 : !! OUTPUT
96 : !! cg2(2, npw2, nspinor, ndat)= symmetrized wavefunctions.
97 : !! work(2, work_ngfft(4), work_ngfft(5), work_ngfft(6))) = workspace array. See comments in INPUTS section.
98 : !!
99 : !! NOTES
100 : !! Inspired to wfconv.
101 : !!
102 : !! SOURCE
103 :
104 481355 : subroutine cgtk_rotate(cryst, kpt1, isym, itimrev, g0, nspinor, ndat, &
105 481355 : npw1, kg1, npw2, kg2, istwf1, istwf2, cg1, cg2, work_ngfft, work)
106 :
107 : !Arguments ------------------------------------
108 : !scalars
109 : integer,intent(in) :: isym, itimrev, nspinor, ndat, npw1, npw2, istwf1, istwf2
110 : type(crystal_t),intent(in) :: cryst
111 : !arrays
112 : integer,intent(in) :: g0(3), kg1(3,npw1), kg2(3,npw2), work_ngfft(18)
113 : real(dp),intent(in) :: kpt1(3), cg1(2,npw1,nspinor,ndat)
114 : real(dp),intent(out) :: cg2(2,npw2,nspinor,ndat)
115 : real(dp),intent(out) :: work(2,work_ngfft(4),work_ngfft(5),work_ngfft(6)) !*ndat) for threads?
116 :
117 : !Local variables ------------------------------
118 : !scalars
119 : integer :: n1,n2,n3,n4,n5,n6,ipw,idat,isp
120 : real(dp) :: arg,ar,ai,bi,br,spinrots,spinrotx,spinroty,spinrotz
121 : logical :: have_phase
122 : !arrays
123 : integer,parameter :: atindx(1) = 1
124 : integer :: symrec(3,3), symrel(3,3)
125 : real(dp) :: phktnons(2,1), tau(3), spinrot(4), tsec(2)
126 481355 : real(dp),allocatable :: phase1d(:,:), phase3d(:,:), wavef1(:,:)
127 : !************************************************************************
128 :
129 : ! Keep track of total time spent.
130 481355 : call timab(1780, 1, tsec)
131 :
132 481355 : ABI_CHECK_IRANGE(itimrev, 0, 1, "itimrev should be in [0, 1]")
133 :
134 481355 : n1 = work_ngfft(1); n2 = work_ngfft(2); n3 = work_ngfft(3)
135 481355 : n4 = work_ngfft(4); n5 = work_ngfft(5); n6 = work_ngfft(6)
136 :
137 6257615 : symrel = cryst%symrel(:,:,isym)
138 481355 : call mati3inv(symrel, symrec) ! symrec = symrel^{-1t}
139 1925420 : tau = cryst%tnons(:,isym)
140 1925420 : have_phase = sum(tau ** 2) > tol8
141 :
142 : ! Compute rotation in spinor space
143 481355 : if (nspinor == 2) call getspinrot(cryst%rprimd, spinrot, symrel)
144 963923 : if (itimrev == 1) symrec = -symrec
145 :
146 : ! Need to compute phase factors associated with nonsymmorphic translations?
147 481355 : if (have_phase) then
148 : ! Although the routine getph is originally written for atomic phase factors, it does precisely what we want
149 40353 : ABI_MALLOC(phase1d, (2, (2*n1+1)+(2*n2+1)+(2*n3+1)))
150 13451 : call getph(atindx, 1, n1, n2, n3, phase1d, tau)
151 :
152 13451 : arg = two_pi * (kpt1(1)*tau(1) + kpt1(2)*tau(2) + kpt1(3)*tau(3))
153 13451 : phktnons(1, 1) = cos(arg)
154 13451 : phktnons(2, 1) = sin(arg)
155 : ! Convert 1D phase factors to 3D phase factors exp(i 2 pi (k1 + G).tnons )
156 40353 : ABI_MALLOC(phase3d, (2, npw1))
157 13451 : call ph1d3d(1, 1, kg1, 1, 1, npw1, n1, n2, n3, phktnons, phase1d, phase3d)
158 13451 : ABI_FREE(phase1d)
159 : end if
160 :
161 1444065 : ABI_MALLOC(wavef1, (2, npw1))
162 :
163 1553713 : do idat=1,ndat
164 2153726 : do isp=1,nspinor
165 1358679450 : wavef1 = cg1(:,:,isp,idat)
166 :
167 1081368 : if (have_phase) then
168 : ! Multiply by phase factors due to nonsymmorphic translations.
169 25635223 : do ipw=1,npw1
170 25573785 : ar = phase3d(1,ipw) * wavef1(1,ipw) - phase3d(2,ipw) * wavef1(2,ipw)
171 25573785 : ai = phase3d(2,ipw) * wavef1(1,ipw) + phase3d(1,ipw) * wavef1(2,ipw)
172 25573785 : wavef1(1, ipw) = ar
173 25635223 : wavef1(2, ipw) = ai
174 : end do
175 : end if
176 :
177 : ! Take into account time-reversal symmetry for SCALAR wavefunctions, if needed.
178 18960484 : if (itimrev == 1 .and. nspinor == 1) wavef1(2, :npw1) = -wavef1(2, :npw1)
179 :
180 : ! Insert wavef1 in work array.
181 1081368 : call sphere(wavef1,ndat1,npw1,work,n1,n2,n3,n4,n5,n6,kg1,istwf1,to_box,me_g0,no_shift,identity_3d,one)
182 :
183 : ! Apply rotation + g0 and extract data on the kg2 sphere: cg2(g) = work(S(g + g0))
184 2153726 : call sphere(cg2(:,:,isp,idat),ndat1,npw2,work,n1,n2,n3,n4,n5,n6,kg2,istwf2,to_sph,me_g0,g0,symrec,one)
185 : end do ! isp
186 :
187 1553713 : if (nspinor == 2) then
188 9010 : if (itimrev == 1) then
189 : ! Take care of time-reversal symmetry, if needed
190 : ! 1) Exchange spin-up and spin-down.
191 : ! 2) Make complex conjugate of one component, and change sign of other component
192 9384 : do ipw=1,npw2
193 : ! Here, change sign of real part
194 9216 : ar = -cg2(1,ipw,1,idat)
195 9216 : ai = cg2(2,ipw,1,idat)
196 : ! Here, change sign of imaginary part
197 9216 : cg2(1,ipw,1,idat) = cg2(1,ipw,2,idat)
198 9216 : cg2(2,ipw,1,idat) = -cg2(2,ipw,2,idat)
199 9216 : cg2(1,ipw,2,idat) = ar
200 9384 : cg2(2,ipw,2,idat) = ai
201 : end do
202 : end if ! itimrev==1
203 :
204 : ! Rotation in spinor space (see also wfconv)
205 9010 : spinrots = spinrot(1); spinrotx = spinrot(2); spinroty = spinrot(3); spinrotz = spinrot(4)
206 2290124 : do ipw=1,npw2
207 2281114 : ar = cg2(1,ipw,1,idat)
208 2281114 : ai = cg2(2,ipw,1,idat)
209 2281114 : br = cg2(1,ipw,2,idat)
210 2281114 : bi = cg2(2,ipw,2,idat)
211 2281114 : cg2(1,ipw,1,idat) = spinrots*ar - spinrotz*ai + spinroty*br - spinrotx*bi
212 2281114 : cg2(2,ipw,1,idat) = spinrots*ai + spinrotz*ar + spinroty*bi + spinrotx*br
213 2281114 : cg2(1,ipw,2,idat) = -spinroty*ar - spinrotx*ai + spinrots*br + spinrotz*bi
214 2290124 : cg2(2,ipw,2,idat) = -spinroty*ai + spinrotx*ar + spinrots*bi - spinrotz*br
215 : end do
216 : end if
217 : end do ! idat
218 :
219 481355 : ABI_FREE(wavef1)
220 481355 : ABI_SFREE(phase3d)
221 :
222 481355 : call timab(1780, 2, tsec)
223 :
224 481355 : end subroutine cgtk_rotate
225 : !!***
226 :
227 : !!****f* ABINIT/cgtk_rotate_symrec
228 : !! NAME
229 : !! cgtk_rotate_symrec
230 : !!
231 : !! FUNCTION
232 : !! Reconstruct wavefunction cg2 in the BZ from the symmetrical image cg1 by applying a symmetry operation.
233 : !! Note that there are two possible conventions for mapping k-points:
234 : !!
235 : !! 1) k2 = T symrel(:,:, isym)^t k1 + g0 (note transpose of symrel)
236 : !!
237 : !! 2) k2 = T symrec(:,:, isym) k1 + g0
238 : !!
239 : !! where T is for time-reversal (itimrev)
240 : !!
241 : !! This routine assumes the SECOND convention.
242 : !!
243 : !! For scalar wavefunctions, we have (with S being a symrec operation)
244 : !!
245 : !! 1) u_{Sk}(g) = e^{-i(Sk + g).tau)} u_k(S^{-1} g) if g0 = 0 and no TR
246 : !!
247 : !! 2) u_{-k}(g) = u_{k}(-g)^* for TR
248 : !!
249 : !! 3) u_{k+g0}(g) = u_{k}(g+g0) if g0 != 0
250 : !!
251 :
252 0 : subroutine cgtk_rotate_symrec(cryst, kpt1, isym, itimrev, g0, nspinor, ndat, &
253 0 : npw1, kg1, npw2, kg2, istwf1, istwf2, cg1, cg2, work_ngfft, work)
254 :
255 : !Arguments ------------------------------------
256 : !scalars
257 : integer,intent(in) :: isym, itimrev, nspinor, ndat, npw1, npw2, istwf1, istwf2
258 : type(crystal_t),intent(in) :: cryst
259 : !arrays
260 : integer,intent(in) :: g0(3), kg1(3,npw1), kg2(3,npw2), work_ngfft(18)
261 : real(dp),intent(in) :: kpt1(3), cg1(2,npw1,nspinor,ndat)
262 : real(dp),intent(out) :: cg2(2,npw2,nspinor,ndat)
263 : real(dp),intent(out) :: work(2,work_ngfft(4),work_ngfft(5),work_ngfft(6)) !*ndat) for threads?
264 :
265 : !Local variables ------------------------------
266 : !scalars
267 : integer :: n1,n2,n3,n4,n5,n6,ipw,idat,isp
268 : real(dp) :: arg,ar,ai,bi,br,spinrots,spinrotx,spinroty,spinrotz
269 : logical :: have_phase
270 : !arrays
271 : integer,parameter :: atindx(1) = 1
272 : integer :: symrec(3,3), symrec_inv(3,3), symrel(3,3), symrel_inv(3,3)
273 : real(dp) :: phktnons(2,1), tau(3), spinrot(4), tsec(2), kpt2(3)
274 0 : real(dp),allocatable :: phase1d(:,:), phase3d(:,:), wavef1(:,:)
275 : !************************************************************************
276 :
277 : ! Keep track of total time spent.
278 0 : call timab(1780, 1, tsec)
279 :
280 0 : ABI_CHECK_IRANGE(itimrev, 0, 1, "itimrev should be in [0, 1]")
281 :
282 0 : n1 = work_ngfft(1); n2 = work_ngfft(2); n3 = work_ngfft(3)
283 0 : n4 = work_ngfft(4); n5 = work_ngfft(5); n6 = work_ngfft(6)
284 :
285 0 : symrec = cryst%symrec(:,:,isym)
286 0 : symrel = cryst%symrel(:,:,isym) ! symrel = symrec^{-1t}
287 : symrec_inv = transpose(symrel)
288 : symrel_inv = transpose(symrec)
289 0 : tau = cryst%tnons(:,isym)
290 0 : have_phase = sum(tau ** 2) > tol8
291 0 : kpt2 = (merge(1, -1, itimrev == 0) * matmul(symrec, kpt1)) + g0
292 :
293 : ! Compute rotation in spinor space
294 0 : if (nspinor == 2) call getspinrot(cryst%rprimd, spinrot, symrel)
295 0 : if (itimrev == 1) symrec = -symrec
296 :
297 : ! Need to compute phase factors associated with nonsymmorphic translations?
298 0 : if (have_phase) then
299 :
300 : ! Although the routine getph is originally written for atomic phase factors, it does precisely what we want
301 0 : ABI_MALLOC(phase1d, (2, (2*n1+1)+(2*n2+1)+(2*n3+1)))
302 0 : call getph(atindx, 1, n1, n2, n3, phase1d, tau)
303 :
304 : ! Note kpt2 instead of kpt1 (difference wrt cgtk_rotate).
305 0 : arg = two_pi * (kpt2(1)*tau(1) + kpt2(2)*tau(2) + kpt2(3)*tau(3))
306 0 : phktnons(1, 1) = cos(arg)
307 0 : phktnons(2, 1) = sin(arg)
308 : ! Convert 1D phase factors to 3D phase factors exp(i 2 pi (k2+G).tnons )
309 0 : ABI_MALLOC(phase3d, (2, npw1))
310 0 : call ph1d3d(1, 1, kg1, 1, 1, npw1, n1, n2, n3, phktnons, phase1d, phase3d)
311 0 : ABI_FREE(phase1d)
312 : end if
313 :
314 0 : ABI_MALLOC(wavef1, (2, npw1))
315 :
316 0 : do idat=1,ndat
317 0 : do isp=1,nspinor
318 0 : wavef1 = cg1(:,:,isp,idat)
319 :
320 0 : if (have_phase) then
321 : ! Multiply by phase factors due to nonsymmorphic translations.
322 : ! Here take complex conjugate of phased3d (note the difference wrt cgtk_rotate).
323 0 : do ipw=1,npw1
324 0 : ar = phase3d(1,ipw) * wavef1(1,ipw) + phase3d(2,ipw) * wavef1(2,ipw)
325 0 : ai = -phase3d(2,ipw) * wavef1(1,ipw) + phase3d(1,ipw) * wavef1(2,ipw)
326 0 : wavef1(1, ipw) = ar
327 0 : wavef1(2, ipw) = ai
328 : end do
329 : end if
330 :
331 : ! Take into account time-reversal symmetry for SCALAR wavefunctions, if needed.
332 0 : if (itimrev == 1 .and. nspinor == 1) wavef1(2, :npw1) = -wavef1(2, :npw1)
333 :
334 : ! Insert wavef1 in work array.
335 0 : call sphere(wavef1,ndat1,npw1,work,n1,n2,n3,n4,n5,n6,kg1,istwf1,to_box,me_g0,no_shift,identity_3d,one)
336 :
337 0 : ABI_ERROR("sphere is not yet compatible with the symrec convention!!")
338 : ! Apply rotation + g0 and extract data on the kg2 sphere: cg2(g) = work(S(g + g0))
339 0 : call sphere(cg2(:,:,isp,idat),ndat1,npw2,work,n1,n2,n3,n4,n5,n6,kg2,istwf2,to_sph,me_g0,g0,symrec,one)
340 : end do ! isp
341 :
342 0 : if (nspinor == 2) then
343 0 : ABI_ERROR("nspinor 2 in cgtk_rotate_symrec is not coded!")
344 0 : if (itimrev == 1) then
345 : ! Take care of time-reversal symmetry, if needed
346 : ! 1) Exchange spin-up and spin-down.
347 : ! 2) Make complex conjugate of one component, and change sign of other component
348 0 : do ipw=1,npw2
349 : ! Here, change sign of real part
350 0 : ar = -cg2(1,ipw,1,idat)
351 0 : ai = cg2(2,ipw,1,idat)
352 : ! Here, change sign of imaginary part
353 0 : cg2(1,ipw,1,idat) = cg2(1,ipw,2,idat)
354 0 : cg2(2,ipw,1,idat) = -cg2(2,ipw,2,idat)
355 0 : cg2(1,ipw,2,idat) = ar
356 0 : cg2(2,ipw,2,idat) = ai
357 : end do
358 : end if ! itimrev==1
359 :
360 : ! Rotation in spinor space (see also wfconv)
361 0 : spinrots = spinrot(1); spinrotx = spinrot(2); spinroty = spinrot(3); spinrotz = spinrot(4)
362 0 : do ipw=1,npw2
363 0 : ar = cg2(1,ipw,1,idat)
364 0 : ai = cg2(2,ipw,1,idat)
365 0 : br = cg2(1,ipw,2,idat)
366 0 : bi = cg2(2,ipw,2,idat)
367 0 : cg2(1,ipw,1,idat) = spinrots*ar - spinrotz*ai + spinroty*br - spinrotx*bi
368 0 : cg2(2,ipw,1,idat) = spinrots*ai + spinrotz*ar + spinroty*bi + spinrotx*br
369 0 : cg2(1,ipw,2,idat) = -spinroty*ar - spinrotx*ai + spinrots*br + spinrotz*bi
370 0 : cg2(2,ipw,2,idat) = -spinroty*ai + spinrotx*ar + spinrots*bi - spinrotz*br
371 : end do
372 : end if
373 : end do ! idat
374 :
375 0 : ABI_FREE(wavef1)
376 0 : ABI_SFREE(phase3d)
377 :
378 0 : call timab(1780, 2, tsec)
379 :
380 0 : end subroutine cgtk_rotate_symrec
381 : !!***
382 :
383 : !!****f* ABINIT/cgtk_change_gsphere
384 : !! NAME
385 : !! cgtk_change_gsphere
386 : !!
387 : !! FUNCTION
388 : !! Transfer the G components of ndat wavefunctions from one sphere to another one.
389 : !! Can also be used to change the value of istwfk e.g. 2 --> 1
390 : !!
391 : !! INPUTS
392 : !! ndat = Number of wavefunctions to transform.
393 : !! npw1, npw2 = Number of plane-waves in the (input, output) G-sphere
394 : !! istwf1, istwf2 = Storage mode of (input, output) wavefunctions.
395 : !! kg1(3,npw1), kg2(3,npw2) = Input/Output G-sphere
396 : !! cg1(2,npw1,ndat) = Input wavefunctions on kg1 sphere with istwf1 mode.
397 : !! work_ngfft(18)=Specify work dimensions. Must be large enough to accommodate kg1 and kg2
398 : !!
399 : !! OUTPUT
400 : !! cg2(2,npw2,ndat) = Output wavefunctions on kg2 sphere with istwf2 mode.
401 : !! work(2,work_ngfft(4),work_ngfft(5),work_ngfft(6)) = Workspace array
402 : !!
403 : !! SOURCE
404 :
405 16946 : subroutine cgtk_change_gsphere(ndat, npw1, istwf1, kg1, cg1, npw2, istwf2, kg2, cg2, work_ngfft, work, &
406 : shiftg1) ! optional
407 :
408 : !Arguments ------------------------------------
409 : !scalars
410 : integer,intent(in) :: ndat,npw1,npw2,istwf1,istwf2
411 : !arrays
412 : integer,intent(in) :: kg1(3,npw1),kg2(3,npw2)
413 : integer,intent(in) :: work_ngfft(18)
414 : real(dp),intent(inout) :: cg1(2,npw1,ndat) ! TODO: Should be intent(in) but need to change sphere
415 : real(dp),intent(out) :: cg2(2,npw2,ndat)
416 : real(dp),intent(out) :: work(2,work_ngfft(4),work_ngfft(5),work_ngfft(6))
417 : integer,optional,intent(in) :: shiftg1(3)
418 :
419 : !Local variables ------------------------------
420 : integer :: n1,n2,n3,n4,n5,n6,idat
421 : integer :: shiftg1__(3)
422 : !************************************************************************
423 :
424 16946 : n1 = work_ngfft(1); n2 = work_ngfft(2); n3 = work_ngfft(3)
425 16946 : n4 = work_ngfft(4); n5 = work_ngfft(5); n6 = work_ngfft(6)
426 :
427 16946 : shiftg1__ = no_shift; if (present(shiftg1)) shiftg1__ = shiftg1
428 :
429 : !print *, "npw1", npw1, "npw2", npw2
430 :
431 35286 : do idat=1,ndat
432 : ! Insert cg1 in work array taking into account istwf1 (intent in)
433 : ! Note: shiftg is only used by sphere when iflag=-1 (extraction).
434 18340 : call sphere(cg1(:,:,idat),1,npw1,work,n1,n2,n3,n4,n5,n6,kg1,istwf1,to_box,me_g0,no_shift,identity_3d,one)
435 :
436 : ! Extract cg2 from work array taking into account istwf2
437 35286 : call sphere(cg2(:,:,idat),1,npw2,work,n1,n2,n3,n4,n5,n6,kg2,istwf2,to_sph,me_g0,shiftg1__,identity_3d,one)
438 : end do
439 :
440 16946 : end subroutine cgtk_change_gsphere
441 : !!***
442 :
443 : !!****f* ABINIT/cgtk_fixphase
444 : !! NAME
445 : !! cgtk_fixphase
446 : !!
447 : !! FUNCTION
448 : !! Fix phase of all bands. Keep normalization but maximize real part
449 : !! (minimize imag part). Also fix the sign of real part
450 : !! by setting the first non-zero element to be positive.
451 : !! See also fxphas_seq in m_cgtools
452 : !!
453 : !! INPUTS
454 : !! cg(2,mcg)= contains the wavefunction |c> coefficients.
455 : !! gsc(2,mgsc)= if useoverlap==1, contains the S|c> coefficients, where S is an overlap matrix.
456 : !! icg=shift to be applied on the location of data in the array cg
457 : !! igsc=shift to be applied on the location of data in the array gsc
458 : !! istwfk=input option parameter that describes the storage of wfs (set to 1 if usual complex vectors)
459 : !! mcg=size of second dimension of cg
460 : !! mgsc=size of second dimension of gsc
461 : !! mpi_enreg=information about MPI parallelization
462 : !! nband_k=number of bands
463 : !! npw_k=number of planewaves
464 : !! useoverlap=describe the overlap of wavefunctions:
465 : !! 0: no overlap (S=Identi0,ty_matrix)
466 : !! 1: wavefunctions are overlapping
467 : !!
468 : !! OUTPUT
469 : !! cg(2,mcg)=same array with altered phase.
470 : !! gsc(2,mgsc)= same array with altered phase.
471 : !!
472 : !! SOURCE
473 :
474 453586 : subroutine cgtk_fixphase(cg, gsc, icg, igsc, istwfk, mcg, mgsc, mpi_enreg, nband_k, npw_k, useoverlap, cprj, nspinor)
475 :
476 : !Arguments ------------------------------------
477 : !scalars
478 : integer,intent(in) :: icg,igsc,istwfk,mcg,mgsc,nband_k,npw_k,useoverlap
479 : type(MPI_type),intent(in) :: mpi_enreg
480 : !arrays
481 : real(dp),intent(inout) :: cg(2,mcg),gsc(2,mgsc*useoverlap)
482 : type(pawcprj_type),intent(inout),optional,target :: cprj(:,:)
483 : integer,intent(in),optional :: nspinor
484 :
485 : !Local variables-------------------------------
486 : !scalars
487 : logical :: do_cprj
488 : integer :: iband,ierr,ii,indx,ncprj
489 : real(dp) :: cim,cre,gscim,gscre,quotient,root1,root2,saa,sab,sbb,theta,thppi,xx,yy
490 : character(len=500) :: msg
491 : !arrays
492 907172 : real(dp) :: buffer2(nband_k,2),buffer3(nband_k,3),tsec(2)
493 453586 : real(dp),allocatable :: cimb(:),creb(:),saab(:),sabb(:),sbbb(:) !,sarr(:,:)
494 : ! *************************************************************************
495 :
496 453586 : do_cprj=.false.
497 453586 : if (present(cprj)) then
498 0 : do_cprj=.true.
499 0 : ncprj = size(cprj,2)
500 0 : if (ncprj/=nband_k*nspinor) then
501 0 : ABI_ERROR('bad size for cprj')
502 : end if
503 : end if
504 :
505 : !The general case, where a complex phase indeterminacy is present
506 453586 : if(istwfk==1)then
507 :
508 1272225 : ABI_MALLOC(cimb,(nband_k))
509 848150 : ABI_MALLOC(creb,(nband_k))
510 848150 : ABI_MALLOC(saab,(nband_k))
511 848150 : ABI_MALLOC(sabb,(nband_k))
512 848150 : ABI_MALLOC(sbbb,(nband_k))
513 6419201 : cimb(:)=zero ; creb(:)=zero
514 :
515 : ! Loop over bands
516 : ! TODO: MG store saa arrays in sarr(3,nband_k) to reduce false sharing.
517 : !$OMP PARALLEL DO DEFAULT(PRIVATE) SHARED(nband_k,icg,npw_k,cg,saab,sbbb,sabb)
518 3421638 : do iband=1,nband_k
519 2997563 : indx=icg+(iband-1)*npw_k
520 :
521 : ! Compute several sums over Re, Im parts of c
522 2997563 : saa=zero; sbb=zero; sab=zero
523 913877405 : do ii=1+indx,npw_k+indx
524 910879842 : saa=saa+cg(1,ii)*cg(1,ii)
525 910879842 : sbb=sbb+cg(2,ii)*cg(2,ii)
526 913877405 : sab=sab+cg(1,ii)*cg(2,ii)
527 : end do
528 2997563 : saab(iband)=saa
529 2997563 : sbbb(iband)=sbb
530 3421638 : sabb(iband)=sab
531 : end do
532 :
533 : ! XG030513 : MPIWF : should transmit saab,sbbb,sabb from the procs
534 : ! of the WF group to the master processor of the WF group
535 424075 : if (mpi_enreg%paral_kgb == 1) then
536 0 : buffer3(:,1)=saab(:)
537 0 : buffer3(:,2)=sbbb(:)
538 0 : buffer3(:,3)=sabb(:)
539 0 : call timab(48,1,tsec)
540 0 : call xmpi_sum(buffer3,mpi_enreg%comm_fft,ierr)
541 0 : if (mpi_enreg%paral_spinor==1) then
542 0 : call xmpi_sum(buffer3,mpi_enreg%comm_spinor,ierr)
543 : end if
544 0 : call timab(48,2,tsec)
545 0 : saab(:)=buffer3(:,1)
546 0 : sbbb(:)=buffer3(:,2)
547 0 : sabb(:)=buffer3(:,3)
548 : end if
549 :
550 : ! XG030513 : MPIWF this loop should only be executed by the master of the WF group
551 :
552 424075 : if (mpi_enreg%paral_kgb==0.or.mpi_enreg%me_fft==0) then
553 3421638 : do iband=1,nband_k
554 2997563 : indx=icg+(iband-1)*npw_k
555 :
556 2997563 : saa=saab(iband)
557 2997563 : sbb=sbbb(iband)
558 2997563 : sab=sabb(iband)
559 :
560 : ! Get phase angle theta
561 2997563 : if (sbb+saa>tol8)then
562 2997563 : if(abs(sbb-saa)>tol8*(sbb+saa) .or. 2*abs(sab)>tol8*(sbb+saa))then
563 2992730 : if (abs(sbb-saa)>tol8*abs(sab)) then
564 2992621 : quotient=sab/(sbb-saa)
565 2992621 : theta=0.5_dp*atan(2.0_dp*quotient)
566 : else
567 : ! Taylor expansion of the atan in terms of inverse of its argument. Correct up to 1/x2, included.
568 109 : theta=0.25_dp*(pi-(sbb-saa)/sab)
569 : end if
570 : ! Check roots to get theta for max Re part
571 2992730 : root1=cos(theta)**2*saa+sin(theta)**2*sbb-2.0_dp*cos(theta)*sin(theta)*sab
572 2992730 : thppi=theta+0.5_dp*pi
573 2992730 : root2=cos(thppi)**2*saa+sin(thppi)**2*sbb-2.0_dp*cos(thppi)*sin(thppi)*sab
574 2992730 : if (root2>root1) theta=thppi
575 : else
576 : ! The real part vector and the imaginary part vector are orthogonal, and of same norm. Strong indeterminacy.
577 : ! Will determine the first non-zero coefficient, and fix its phase
578 : ! Hypothesis : there is at least one non-zero element on the master node ...
579 397414 : do ii=1+indx,npw_k+indx
580 397414 : cre=cg(1,ii)
581 397414 : cim=cg(2,ii)
582 397414 : if(cre**2+cim**2>tol8**2*(saa+sbb))then
583 4833 : if(cre**2>tol8**2**cim**2)then
584 3941 : theta=atan(cim/cre)
585 : else
586 : ! Taylor expansion of the atan in terms of inverse of its argument. Correct up to 1/x2, included.
587 892 : theta=pi/2-cre/cim
588 : end if
589 : exit
590 : end if
591 : end do
592 : end if
593 : else
594 : write(msg,'(a,i0,5a)')&
595 0 : & 'The eigenvector with band ',iband,' has zero norm.',ch10,&
596 0 : & 'This usually happens when the number of bands (nband) is comparable to the number of planewaves (mpw)',ch10,&
597 0 : & 'Action: Check the parameters of the calculation. If nband ~ mpw, then decrease nband or, alternatively, increase ecut'
598 0 : ABI_ERROR(msg)
599 : end if
600 :
601 2997563 : xx=cos(theta)
602 2997563 : yy=sin(theta)
603 :
604 : ! Here, set the first non-zero element to be positive
605 : ! Comment the next nine lines to recover the behaviour of pre v3.1.3g
606 : ! Hypothesis : there is at least one non-zero element on the master node ...
607 7088078 : do ii=1+indx,npw_k+indx
608 7088078 : cre=cg(1,ii)
609 7088078 : cim=cg(2,ii)
610 7088078 : cre=xx*cre-yy*cim
611 7088078 : if(abs(cre)>tol8)exit
612 : end do
613 2997563 : if(cre<zero)then
614 411836 : xx=-xx ; yy=-yy
615 : end if
616 :
617 2997563 : creb(iband)=xx
618 3421638 : cimb(iband)=yy
619 :
620 : end do
621 : end if
622 :
623 : ! XG030513 : MPIWF : should transmit creb(:),cimb(:) of the master
624 : ! processor of the WF group to the others procs of the WF group
625 424075 : if (mpi_enreg%paral_kgb == 1) then
626 0 : call timab(48,1,tsec)
627 0 : buffer2(:,1)=creb(:)
628 0 : buffer2(:,2)=cimb(:)
629 0 : call xmpi_sum(buffer2,mpi_enreg%comm_fft,ierr)
630 0 : if (mpi_enreg%paral_spinor==1) then
631 0 : call xmpi_sum(buffer2,mpi_enreg%comm_spinor,ierr)
632 : end if
633 0 : call timab(48,2,tsec)
634 0 : creb(:)=buffer2(:,1)
635 0 : cimb(:)=buffer2(:,2)
636 : end if
637 :
638 : ! MG TODO: Scaling can be done with zscal
639 : !$OMP PARALLEL DO PRIVATE(indx,xx,yy,cre,cim,gscre,gscim)
640 3421638 : do iband=1,nband_k
641 2997563 : indx=icg+(iband-1)*npw_k
642 :
643 2997563 : xx=creb(iband)
644 2997563 : yy=cimb(iband)
645 : ! Alter phase of array |cg>
646 913877405 : do ii=1+indx,npw_k+indx
647 910879842 : cre=cg(1,ii)
648 910879842 : cim=cg(2,ii)
649 910879842 : cg(1,ii)=xx*cre-yy*cim
650 913877405 : cg(2,ii)=xx*cim+yy*cre
651 : end do
652 2997563 : if (do_cprj) call pawcprj_zaxpby((/zero,zero/),(/xx,yy/),cprj(:,nspinor*(iband-1)+1:nspinor*iband),&
653 0 : & cprj(:,nspinor*(iband-1)+1:nspinor*iband))
654 :
655 : ! Alter phase of array S|cg>
656 3421638 : if (useoverlap==1) then
657 1172917 : indx=igsc+(iband-1)*npw_k
658 370101406 : do ii=1+indx,npw_k+indx
659 368928489 : gscre=gsc(1,ii)
660 368928489 : gscim=gsc(2,ii)
661 368928489 : gsc(1,ii)=xx*gscre-yy*gscim
662 370101406 : gsc(2,ii)=xx*gscim+yy*gscre
663 : end do
664 : end if
665 : end do ! iband
666 :
667 424075 : ABI_FREE(cimb)
668 424075 : ABI_FREE(creb)
669 424075 : ABI_FREE(saab)
670 424075 : ABI_FREE(sabb)
671 424075 : ABI_FREE(sbbb)
672 :
673 : else ! if istwfk/=1. Storages that take into account the time-reversal symmetry : the freedom is only a sign freedom
674 :
675 88533 : ABI_MALLOC(creb,(nband_k))
676 294231 : creb(:)=zero
677 : ! XG030513 : MPIWF : this loop should be done only by the master processor of the WF group
678 :
679 29511 : if (mpi_enreg%paral_kgb==0.or.mpi_enreg%me_fft==0) then
680 :
681 : ! Loop over bands
682 294231 : do iband=1,nband_k
683 :
684 264720 : indx=icg+(iband-1)*npw_k
685 :
686 : ! Here, set the first non-zero real element to be positive
687 10769542 : do ii=1+indx,npw_k+indx
688 10749844 : cre=cg(1,ii)
689 10769542 : if(abs(cre)>tol8)exit
690 : end do
691 294231 : creb(iband)=cre
692 :
693 : end do ! iband
694 :
695 : end if
696 : ! XG030513 : MPIWF : should transmit cre(:) of the master processor of the WF group to the others
697 29511 : if (mpi_enreg%paral_kgb == 1) then
698 0 : call timab(48,1,tsec)
699 0 : call xmpi_sum(creb,mpi_enreg%comm_fft,ierr)
700 0 : if (mpi_enreg%paral_spinor==1) then
701 0 : call xmpi_sum(creb,mpi_enreg%comm_spinor,ierr)
702 : end if
703 0 : call timab(48,2,tsec)
704 : end if
705 :
706 294231 : do iband=1,nband_k
707 264720 : cre=creb(iband)
708 294231 : if(cre<zero)then
709 51269 : indx=icg+(iband-1)*npw_k
710 25275750 : do ii=1+indx,npw_k+indx
711 25224481 : cg(1,ii)=-cg(1,ii)
712 25275750 : cg(2,ii)=-cg(2,ii)
713 : end do
714 51269 : if (do_cprj) call pawcprj_zaxpby((/zero,zero/),(/-one,zero/),cprj(:,iband:iband),cprj(:,iband:iband))
715 51269 : if(useoverlap==1)then
716 14417 : indx=igsc+(iband-1)*npw_k
717 11239482 : do ii=1+indx,npw_k+indx
718 11225065 : gsc(1,ii)=-gsc(1,ii)
719 11239482 : gsc(2,ii)=-gsc(2,ii)
720 : end do
721 : end if
722 : end if
723 : end do ! iband
724 :
725 29511 : ABI_FREE(creb)
726 : end if ! istwfk
727 :
728 453586 : end subroutine cgtk_fixphase
729 : !!***
730 :
731 : end module m_cgtk
732 : !!***
|
