Branch data 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 : : !!***
|
