Branch data Line data Source code
1 : : !!****m* ABINIT/m_spacepar
2 : : !! NAME
3 : : !! m_spacepar
4 : : !!
5 : : !! FUNCTION
6 : : !! Relatively Low-level procedures operating on arrays defined on the FFT box (G- or R- space)
7 : : !! Unlike the procedures in m_cgtools, the routines declared in this module can use mpi_type.
8 : : !!
9 : : !! COPYRIGHT
10 : : !! Copyright (C) 2008-2026 ABINIT group (XG, BA, MT, DRH, DCA, GMR, MJV, JWZ)
11 : : !! This file is distributed under the terms of the
12 : : !! GNU General Public License, see ~abinit/COPYING
13 : : !! or http://www.gnu.org/copyleft/gpl.txt .
14 : : !!
15 : : !! SOURCE
16 : :
17 : : #if defined HAVE_CONFIG_H
18 : : #include "config.h"
19 : : #endif
20 : :
21 : : #include "abi_common.h"
22 : :
23 : : module m_spacepar
24 : :
25 : : use, intrinsic :: iso_c_binding, only: c_loc, c_associated
26 : :
27 : : use defs_basis
28 : : use m_abicore
29 : : use m_errors
30 : : use m_xmpi
31 : : use m_xomp
32 : : use m_sort
33 : :
34 : : use m_time, only : timab
35 : : use defs_abitypes, only : MPI_type
36 : : use m_matrix, only : mati3inv, matr3inv
37 : : use m_symtk, only : sg_multable, symdet, symatm
38 : : use m_geometry, only : metric, normv, symredcart,wedge_basis,wedge_product
39 : : use m_gtermcutoff, only : termcutoff
40 : : use m_mpinfo, only : ptabs_fourdp
41 : : use m_fft, only : zerosym, fourdp
42 : :
43 : : implicit none
44 : :
45 : : private
46 : : !!***
47 : :
48 : : public :: hartre ! Given rho(G), compute Hartree potential (=FFT of rho(G)/pi/(G+q)**2)
49 : : public :: mkunitpawspherepot ! compute effective potential due to PAW sphere of strength 1, useful for testing
50 : : ! paw projector completeness
51 : : public :: make_vectornd ! compute vector potential due to nuclear magnetic dipoles, in real space
52 : : public :: meanvalue_g ! Compute <wf|op|wf> where op is real and diagonal in G-space.
53 : : public :: meanvalue_g_batch ! Batched GPU-accelerated variant of meanvalue_g (istwf_k==1, filter==0, use_ndo==0).
54 : : public :: laplacian ! Compute the laplacian of a function defined in real space
55 : : public :: redgr ! Compute reduced gradients of a real function on the usual unshifted FFT grid.
56 : : public :: hartrestr ! FFT of (rho(G)/pi)*[d(1/G**2)/d(strain) - delta(diagonal strain)*(1/G**2)]
57 : : public :: symrhg ! Symmetrize rhor(r)
58 : : public :: irrzg ! Find the irreducible zone in reciprocal space (used by symrhg)
59 : : public :: setsym ! Set up irreducible zone in G space by direct calculation.
60 : : public :: hartredq ! Compute the q-gradient of the Hartree potential (=FFT of -rho(G)*G_qdir/pi**2/|G|**4 )
61 : :
62 : : ! MG FIXME This routine is deprecated. Now the symmetrization of the DFPT **potentials** is done in the m_dvdb
63 : : public :: rotate_rho
64 : : !!***
65 : :
66 : : contains
67 : : !!***
68 : :
69 : : !!****f* m_spacepar/make_vectornd
70 : : !! NAME
71 : : !! make_vectornd
72 : : !!
73 : : !! FUNCTION
74 : : !! For nuclear dipole moments m, compute vector potential A(r) = \alpha^2(m x (r-R))/|r-R|^3
75 : : !! in r space. This is done by computing A(G) followed by FFT.
76 : : !!
77 : : !! NOTES
78 : : !! This code is copied and modified from m_spacepar/hartre where a very similar loop
79 : : !! over G is done followed by FFT to real space
80 : : !!
81 : : !! INPUTS
82 : : !!
83 : : !! OUTPUT
84 : : !! vectornd(3,nfft)=Vector potential in real space, along Cartesian directions
85 : : !!
86 : : !! SOURCE
87 : :
88 : 110 : subroutine make_vectornd(cplex,gsqcut,izero,mpi_enreg,natom,nfft,ngfft,nspden,nucdipmom,&
89 : 110 : & rprimd,vectornd,xred)
90 : :
91 : : !Arguments ------------------------------------
92 : : !scalars
93 : : integer,intent(in) :: cplex,izero,natom,nfft,nspden
94 : : real(dp),intent(in) :: gsqcut
95 : : type(MPI_type),intent(in) :: mpi_enreg
96 : : !arrays
97 : : integer,intent(in) :: ngfft(18)
98 : : real(dp),intent(in) :: nucdipmom(3,natom),rprimd(3,3),xred(3,natom)
99 : : real(dp),intent(out) :: vectornd(nfft,nspden,3)
100 : :
101 : : !Local variables-------------------------------
102 : : !scalars
103 : : integer,parameter :: im=2,re=1
104 : : integer :: i1,i2,i2_local,i23,i3,iatom,id1,id2,id3,ig,ig1,ig2,ig3,ig1max,ig2max,ig3max
105 : : integer :: ig1min,ig2min,ig3min
106 : : integer :: ii,ii1,ing,me_fft,n1,n2,n3,nd_atom,nd_atom_tot,nproc_fft
107 : : real(dp),parameter :: tolfix=1.000000001e0_dp
108 : : real(dp) :: cutoff,gqgm12,gqg2p3,gqgm23,gqgm13,gs2,gs3,gs,phase,ucvol
109 : : complex(dp) :: prefac,cgr
110 : : !arrays
111 : : integer :: id(3)
112 : 110 : integer,allocatable :: nd_list(:)
113 : 110 : integer, contiguous, pointer :: fftn2_distrib(:),ffti2_local(:)
114 : 110 : integer, contiguous, pointer :: fftn3_distrib(:),ffti3_local(:)
115 : : real(dp) :: gmet(3,3),gprimd(3,3),gqred(3),mcgc(3),rmet(3,3)
116 : : real(dp) :: rgbasis(3,3,3)
117 : 110 : real(dp),allocatable :: gq(:,:),nd_m(:,:),ndvecr(:),work1(:,:),work2(:,:),work3(:,:)
118 : : ! *************************************************************************
119 : :
120 : 110 : call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
121 : :
122 : : ! make list of atoms with nonzero nuclear dipole moments
123 : : ! in typical applications only 0 or 1 atoms have nonzero dipoles. This
124 : : ! code shouldn't even be called if all dipoles are zero.
125 : 110 : nd_atom_tot = 0
126 [ + + ]: 290 : do iatom = 1, natom
127 [ + + + + : 504 : if (any(abs(nucdipmom(:,iatom))>tol8)) then
+ + ]
128 : 116 : nd_atom_tot = nd_atom_tot + 1
129 : : end if
130 : : end do
131 : :
132 : : ! construct the basis vectors of the generalized cross product
133 : : ! real space a, b, c (contained in rprimd)
134 : : ! reciprocal space a*, b*, c* (contained in gprimd)
135 : : ! for m x G will need a x a*, a x b* etc (9 a^b type basis vectors)
136 : 110 : call wedge_basis(gprimd,rprimd,rgbasis)
137 : :
138 : : ! note that nucdipmom is input as vectors in atomic units referenced
139 : : ! to cartesian coordinates
140 [ - + + - : 330 : ABI_MALLOC(nd_list,(nd_atom_tot))
+ - + - -
+ - + ]
141 [ - + + - : 330 : ABI_MALLOC(nd_m,(3,nd_atom_tot))
+ - + - -
+ - + ]
142 : 290 : nd_atom_tot = 0
143 [ + + ]: 290 : do iatom = 1, natom
144 [ + + + + : 504 : if (any(abs(nucdipmom(:,iatom))>tol8)) then
+ + ]
145 : 116 : nd_atom_tot = nd_atom_tot + 1
146 : 116 : nd_list(nd_atom_tot) = iatom
147 : : ! the following expresses the dipole moment components in units of rprimd translations
148 : 116 : nd_m(:,nd_atom_tot) = MATMUL(TRANSPOSE(gprimd),nucdipmom(:,iatom))
149 : : end if
150 : : end do
151 : :
152 : 110 : n1=ngfft(1); n2=ngfft(2); n3=ngfft(3)
153 : 110 : nproc_fft = mpi_enreg%nproc_fft; me_fft = mpi_enreg%me_fft
154 : :
155 : : ! the two_pi in the denominator arises from using G.G=2\pi gmet below
156 : 110 : prefac = -four_pi*j_dpc*FineStructureConstant2/(ucvol*two_pi)
157 : :
158 : : ! Get the distrib associated with this fft_grid
159 : 110 : call ptabs_fourdp(mpi_enreg,n2,n3,fftn2_distrib,ffti2_local,fftn3_distrib,ffti3_local)
160 : :
161 : : ! Initialize a few quantities
162 : 110 : cutoff=gsqcut*tolfix
163 : :
164 : : ! In order to speed the routine, precompute the components of g+q
165 : : ! Also check if the booked space was large enough...
166 [ - + + - : 330 : ABI_MALLOC(gq,(3,max(n1,n2,n3)))
+ - + - -
+ - + ]
167 [ + + ]: 440 : do ii=1,3
168 : 330 : id(ii)=ngfft(ii)/2+2
169 [ + + ]: 15488 : do ing=1,ngfft(ii)
170 : 15048 : ig=ing-(ing/id(ii))*ngfft(ii)-1
171 : 15378 : gq(ii,ing)=ig
172 : : end do
173 : : end do
174 : 110 : ig1max=-1;ig2max=-1;ig3max=-1
175 : 110 : ig1min=n1;ig2min=n2;ig3min=n3
176 : :
177 [ - + + - : 330 : ABI_MALLOC(work1,(2,nfft))
+ - + - -
+ - + ]
178 [ - + + - : 220 : ABI_MALLOC(work2,(2,nfft))
- + - + ]
179 [ - + + - : 220 : ABI_MALLOC(work3,(2,nfft))
- + - + ]
180 [ + + + + : 247158254 : work1=zero; work2=zero; work3=zero
+ + + + +
+ + + ]
181 : 5126 : id1=n1/2+2;id2=n2/2+2;id3=n3/2+2
182 : :
183 : : ! Triple loop on each dimension
184 [ + + ]: 5126 : do i3=1,n3
185 : 5016 : ig3=i3-(i3/id3)*n3-1
186 : : ! Precompute some products that do not depend on i2 and i1
187 : 5016 : gs3=gq(3,i3)*gq(3,i3)*gmet(3,3)
188 : 5016 : gqgm23=gq(3,i3)*gmet(2,3)*2
189 : 5016 : gqgm13=gq(3,i3)*gmet(1,3)*2
190 : :
191 [ + + ]: 348230 : do i2=1,n2
192 : 343104 : ig2=i2-(i2/id2)*n2-1
193 [ + + ]: 348120 : if (fftn2_distrib(i2) == me_fft) then
194 : 323104 : gs2=gs3+ gq(2,i2)*(gq(2,i2)*gmet(2,2)+gqgm23)
195 : 323104 : gqgm12=gq(2,i2)*gmet(1,2)*2
196 : 323104 : gqg2p3=gqgm13+gqgm12
197 : :
198 : 323104 : i2_local = ffti2_local(i2)
199 : 323104 : i23=n1*(i2_local-1 +(n2/nproc_fft)*(i3-1))
200 : : ! Do the test that eliminates the Gamma point outside of the inner loop
201 : 323104 : ii1=1
202 : : !if(i23==0 .and. ig2==0 .and. ig3==0)then
203 : : ! ii1=2
204 : : ! work1(re,1+i23)=zero
205 : : ! work1(im,1+i23)=zero
206 : : !end if
207 : :
208 : : ! Final inner loop on the first dimension (note the lower limit)
209 [ + + ]: 27785120 : do i1=ii1,n1
210 : 27462016 : gs=gs2+ gq(1,i1)*(gq(1,i1)*gmet(1,1)+gqg2p3)
211 : 27462016 : ig1 = i1 - (i1/id1)*n1 -1
212 : 27462016 : ii=i1+i23
213 : :
214 : 27462016 : gqred(1) = gq(1,i1); gqred(2) = gq(2,i2); gqred(3) = gq(3,i3)
215 : :
216 [ + + + + ]: 27785120 : if( (gs .LE. cutoff) .AND. (gs .gt. tol8) )then
217 : :
218 [ + + ]: 27813736 : do iatom = 1, nd_atom_tot
219 : 14464612 : nd_atom = nd_list(iatom)
220 [ + + ]: 57858448 : phase = -two_pi*DOT_PRODUCT(xred(:,nd_atom),gqred(:))
221 : 14464612 : cgr = cmplx(cos(phase),sin(phase))
222 : :
223 : : ! cross product m x G
224 : 14464612 : call wedge_product(mcgc,nd_m(:,iatom),gqred,rgbasis)
225 : :
226 : : ! express mcgc relative to rprimd translations. This is done because
227 : : ! we wish ultimately to apply A.p to |cwavef>; in getghc_nucdip, the
228 : : ! p|cwavef> is done in reduced coordinates so do that here too, because
229 : : ! r.G has no need of the metric if both terms are in reduced coords
230 : 14464612 : mcgc = MATMUL(TRANSPOSE(gprimd),mcgc)
231 : :
232 : 14464612 : work1(re,ii) = work1(re,ii) + real(prefac*cgr*mcgc(1)/gs)
233 : 14464612 : work2(re,ii) = work2(re,ii) + real(prefac*cgr*mcgc(2)/gs)
234 : 14464612 : work3(re,ii) = work3(re,ii) + real(prefac*cgr*mcgc(3)/gs)
235 : :
236 : 14464612 : work1(im,ii) = work1(im,ii) + aimag(prefac*cgr*mcgc(1)/gs)
237 : 14464612 : work2(im,ii) = work2(im,ii) + aimag(prefac*cgr*mcgc(2)/gs)
238 : 27813736 : work3(im,ii) = work3(im,ii) + aimag(prefac*cgr*mcgc(3)/gs)
239 : :
240 : : end do
241 : : else
242 : : ! gs>cutoff
243 : 14112892 : work1(re,ii)=zero
244 : 14112892 : work1(im,ii)=zero
245 : 14112892 : work2(re,ii)=zero
246 : 14112892 : work2(im,ii)=zero
247 : 14112892 : work3(re,ii)=zero
248 : 14112892 : work3(im,ii)=zero
249 : : end if
250 : :
251 : : end do ! End loop on i1
252 : : end if
253 : : end do ! End loop on i2
254 : : end do ! End loop on i3
255 : :
256 [ - + ]: 110 : ABI_FREE(gq)
257 [ - + ]: 110 : ABI_FREE(nd_list)
258 [ - + ]: 110 : ABI_FREE(nd_m)
259 : :
260 [ + + ]: 110 : if ( izero .EQ. 1 ) then
261 : : ! Set contribution of unbalanced components to zero
262 : :
263 : 105 : call zerosym(work1,2,n1,n2,n3,comm_fft=mpi_enreg%comm_fft,distribfft=mpi_enreg%distribfft)
264 : 105 : call zerosym(work2,2,n1,n2,n3,comm_fft=mpi_enreg%comm_fft,distribfft=mpi_enreg%distribfft)
265 : 105 : call zerosym(work3,2,n1,n2,n3,comm_fft=mpi_enreg%comm_fft,distribfft=mpi_enreg%distribfft)
266 : :
267 : : end if
268 : :
269 : : !note nspden effect--the nuclear vector potential conains no electron spin flip operator,
270 : : ! so vectornd(:,2,:) = vectornd(:,1,:) and vectornd(:,3:4,:) = zero
271 [ + + + + : 82635668 : vectornd = zero
+ + ]
272 : : ! Fourier Transform
273 [ - + + - : 330 : ABI_MALLOC(ndvecr,(cplex*nfft))
+ - + - -
+ - + ]
274 [ + + ]: 27462126 : ndvecr=zero
275 : 110 : call fourdp(cplex,work1,ndvecr,1,mpi_enreg,nfft,1,ngfft,0)
276 [ + + ]: 27462126 : vectornd(:,1,1)=ndvecr(:)
277 [ + + + + ]: 83054 : if (nspden .GE. 2) vectornd(:,2,1) = ndvecr(:)
278 [ - + ]: 110 : ABI_FREE(work1)
279 : :
280 [ + + ]: 27462126 : ndvecr=zero
281 : 110 : call fourdp(cplex,work2,ndvecr,1,mpi_enreg,nfft,1,ngfft,0)
282 [ + + ]: 27462126 : vectornd(:,1,2) = ndvecr(:)
283 [ + + + + ]: 83054 : if (nspden .GE. 2) vectornd(:,2,2) = ndvecr(:)
284 [ - + ]: 110 : ABI_FREE(work2)
285 : :
286 [ + + ]: 27462126 : ndvecr=zero
287 : 110 : call fourdp(cplex,work3,ndvecr,1,mpi_enreg,nfft,1,ngfft,0)
288 [ + + ]: 27462126 : vectornd(:,1,3) = ndvecr(:)
289 [ + + + + ]: 83054 : if (nspden .GE. 2) vectornd(:,2,3) = ndvecr(:)
290 [ - + ]: 110 : ABI_FREE(work3)
291 [ - + ]: 110 : ABI_FREE(ndvecr)
292 : :
293 : 110 : end subroutine make_vectornd
294 : : !!***
295 : :
296 : : !!****f* m_spacepar/mkunitpawspherepot
297 : : !! NAME
298 : : !! mkunitpawspherepot
299 : : !!
300 : : !! FUNCTION
301 : : !! Compute "potential" due to a sphere of radius r_paw at one of the ions, of
302 : : !! strength 1. This is done for testing the completeness of the PAW projectors.
303 : : !!
304 : : !! NOTES
305 : : !!
306 : : !! INPUTS
307 : : !!
308 : : !! OUTPUT
309 : : !! vunitpawspherepot(cplex*nfft)=Hartree potential in real space, either REAL or COMPLEX
310 : : !!
311 : : !! SOURCE
312 : :
313 : 0 : subroutine mkunitpawspherepot(cplex,gsqcut,izero,mpi_enreg,natom,nfft,ngfft,&
314 : 0 : & rpaw,rprimd,vunitpawspherepot,xred,&
315 : : the_atom) ! Optional arguments
316 : :
317 : : !Arguments ------------------------------------
318 : : !scalars
319 : : integer,intent(in) :: cplex,izero,natom,nfft
320 : : integer,intent(in),optional :: the_atom
321 : : real(dp),intent(in) :: gsqcut,rpaw
322 : : type(MPI_type),intent(in) :: mpi_enreg
323 : : !arrays
324 : : integer,intent(in) :: ngfft(18)
325 : : real(dp),intent(in) :: rprimd(3,3),xred(3,natom)
326 : : real(dp),intent(out) :: vunitpawspherepot(cplex*nfft)
327 : :
328 : : !Local variables-------------------------------
329 : : !scalars
330 : : integer,parameter :: im=2,re=1
331 : : integer :: i1,i2,i2_local,i23,i3,iatom,id1,id2,id3,ig,ig1,ig2,ig3,ig1max,ig2max,ig3max
332 : : integer :: ig1min,ig2min,ig3min
333 : : integer :: ii,ii1,ing,me_fft,n1,n2,n3,nproc_fft,qeq0,qeq05
334 : : real(dp),parameter :: tolfix=1.000000001e0_dp
335 : : real(dp) :: cutoff,gqgm12,gqg2p3,gqgm23,gqgm13,gs2,gs3,gs,ogg0,ogr,ogrpre,phgr,phpaw,ucvol
336 : : character(len=500) :: message
337 : : !arrays
338 : : integer :: id(3)
339 : 0 : integer, contiguous, pointer :: fftn2_distrib(:),ffti2_local(:)
340 : 0 : integer, contiguous, pointer :: fftn3_distrib(:),ffti3_local(:)
341 : : real(dp) :: gmet(3,3),gprimd(3,3),gqred(3),qpt_(3),rmet(3,3)
342 : 0 : real(dp),allocatable :: gq(:,:),work1(:,:)
343 : : ! *************************************************************************
344 : :
345 : : ! Check that cplex has an allowed value
346 : 0 : if(cplex/=1 .and. cplex/=2)then
347 : : write(message, '(a,i0,a,a)' )&
348 : 0 : 'From the calling routine, cplex=',cplex,ch10, 'but the only value allowed are 1 and 2.'
349 : 0 : ABI_BUG(message)
350 : : end if
351 : :
352 : 0 : call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
353 : :
354 : 0 : ogrpre = four_pi/(ucvol*(two_pi**3))
355 : 0 : ogg0 = four_pi*(rpaw**3)/(three*ucvol)
356 : :
357 : 0 : n1=ngfft(1); n2=ngfft(2); n3=ngfft(3)
358 : 0 : nproc_fft = mpi_enreg%nproc_fft; me_fft = mpi_enreg%me_fft
359 : :
360 : : ! Get the distrib associated with this fft_grid
361 : 0 : call ptabs_fourdp(mpi_enreg,n2,n3,fftn2_distrib,ffti2_local,fftn3_distrib,ffti3_local)
362 : :
363 : : ! Initialize a few quantities
364 : 0 : cutoff=gsqcut*tolfix
365 : : !carrying code over from hartre with minimal changes, in present case qpt always zero
366 : 0 : qpt_=zero
367 : : qeq0=0
368 : : if(qpt_(1)**2+qpt_(2)**2+qpt_(3)**2<1.d-15) qeq0=1
369 : : qeq05=0
370 : : if (qeq0==0) then
371 : : if (abs(abs(qpt_(1))-half)<tol12.or.abs(abs(qpt_(2))-half)<tol12.or.abs(abs(qpt_(3))-half)<tol12) qeq05=1
372 : : end if
373 : :
374 [ # # ]: 0 : if (present(the_atom)) then
375 : 0 : iatom = the_atom
376 : : else
377 : : iatom = 1
378 : : end if
379 : :
380 : : ! If cplex=1 then qpt_ should be 0 0 0
381 : : if (cplex==1.and. qeq0/=1) then
382 : : write(message,'(a,3e12.4,a,a)')&
383 : : 'cplex=1 but qpt=',qpt_,ch10,&
384 : : 'qpt should be 0 0 0.'
385 : : ABI_BUG(message)
386 : : end if
387 : :
388 : : ! If FFT parallelism then qpt should not be 1/2
389 : : if (nproc_fft>1.and.qeq05==1) then
390 : : write(message, '(a,3e12.4,a,a)' )&
391 : : 'FFT parallelism selected but qpt',qpt_,ch10,&
392 : : 'qpt(i) should not be 1/2...'
393 : : ABI_ERROR(message)
394 : : end if
395 : :
396 : : ! In order to speed the routine, precompute the components of g+q
397 : : ! Also check if the booked space was large enough...
398 [ # # # # : 0 : ABI_MALLOC(gq,(3,max(n1,n2,n3)))
# # # # #
# # # ]
399 [ # # ]: 0 : do ii=1,3
400 : 0 : id(ii)=ngfft(ii)/2+2
401 [ # # ]: 0 : do ing=1,ngfft(ii)
402 : 0 : ig=ing-(ing/id(ii))*ngfft(ii)-1
403 : 0 : gq(ii,ing)=ig+qpt_(ii)
404 : : end do
405 : : end do
406 : 0 : ig1max=-1;ig2max=-1;ig3max=-1
407 : 0 : ig1min=n1;ig2min=n2;ig3min=n3
408 : :
409 [ # # # # : 0 : ABI_MALLOC(work1,(2,nfft))
# # # # #
# # # ]
410 [ # # # # ]: 0 : work1=zero
411 : 0 : id1=n1/2+2;id2=n2/2+2;id3=n3/2+2
412 : :
413 : : ! Triple loop on each dimension
414 [ # # ]: 0 : do i3=1,n3
415 : 0 : ig3=i3-(i3/id3)*n3-1
416 : : ! Precompute some products that do not depend on i2 and i1
417 : 0 : gs3=gq(3,i3)*gq(3,i3)*gmet(3,3)
418 : 0 : gqgm23=gq(3,i3)*gmet(2,3)*2
419 : 0 : gqgm13=gq(3,i3)*gmet(1,3)*2
420 : :
421 [ # # ]: 0 : do i2=1,n2
422 : 0 : ig2=i2-(i2/id2)*n2-1
423 [ # # ]: 0 : if (fftn2_distrib(i2) == me_fft) then
424 : 0 : gs2=gs3+ gq(2,i2)*(gq(2,i2)*gmet(2,2)+gqgm23)
425 : 0 : gqgm12=gq(2,i2)*gmet(1,2)*2
426 : 0 : gqg2p3=gqgm13+gqgm12
427 : :
428 : 0 : i2_local = ffti2_local(i2)
429 : 0 : i23=n1*(i2_local-1 +(n2/nproc_fft)*(i3-1))
430 : : ! Do the test that eliminates the Gamma point outside of the inner loop
431 : 0 : ii1=1
432 [ # # # # ]: 0 : if(i23==0 .and. qeq0==1 .and. ig2==0 .and. ig3==0)then
433 : 0 : ii1=2
434 : 0 : work1(re,1+i23)=ogg0
435 : 0 : work1(im,1+i23)=zero
436 : : end if
437 : :
438 : : ! Final inner loop on the first dimension (note the lower limit)
439 [ # # ]: 0 : do i1=ii1,n1
440 : 0 : gs=gs2+ gq(1,i1)*(gq(1,i1)*gmet(1,1)+gqg2p3)
441 : 0 : ii=i1+i23
442 : :
443 : 0 : gqred(1) = gq(1,i1); gqred(2) = gq(2,i2); gqred(3) = gq(3,i3)
444 : :
445 [ # # ]: 0 : if(gs<=cutoff)then
446 : :
447 : :
448 : : ! Identify min/max indexes (to cancel unbalanced contributions later)
449 : : ! Count (q+g)-vectors with similar norm
450 : : if ((qeq05==1).and.(izero==1)) then
451 : : ig1=i1-(i1/id1)*n1-1
452 : : ig1max=max(ig1max,ig1); ig1min=min(ig1min,ig1)
453 : : ig2max=max(ig2max,ig2); ig2min=min(ig2min,ig2)
454 : : ig3max=max(ig3max,ig3); ig3min=min(ig3min,ig3)
455 : : end if
456 : :
457 : 0 : phpaw=two_pi*rpaw*gs
458 : 0 : ogr=ogrpre*(sin(phpaw)-phpaw*cos(phpaw))/(gs**3)
459 : :
460 [ # # ]: 0 : phgr = -two_pi*DOT_PRODUCT(xred(:,iatom),gqred(:))
461 : :
462 : 0 : work1(re,ii)=cos(phgr)*ogr
463 : 0 : work1(im,ii)=sin(phgr)*ogr
464 : : else
465 : : ! gs>cutoff
466 : 0 : work1(re,ii)=zero
467 : 0 : work1(im,ii)=zero
468 : : end if
469 : :
470 : : end do ! End loop on i1
471 : : end if
472 : : end do ! End loop on i2
473 : : end do ! End loop on i3
474 : :
475 [ # # ]: 0 : ABI_FREE(gq)
476 : :
477 [ # # ]: 0 : if (izero==1) then
478 : : ! Set contribution of unbalanced components to zero
479 : :
480 : : if (qeq0==1) then !q=0
481 : 0 : call zerosym(work1,2,n1,n2,n3,comm_fft=mpi_enreg%comm_fft,distribfft=mpi_enreg%distribfft)
482 : :
483 : : else if (qeq05==1) then
484 : : !q=1/2; this doesn't work in parallel
485 : : ig1=-1;if (mod(n1,2)==0) ig1=1+n1/2
486 : : ig2=-1;if (mod(n2,2)==0) ig2=1+n2/2
487 : : ig3=-1;if (mod(n3,2)==0) ig3=1+n3/2
488 : : if (abs(abs(qpt_(1))-half)<tol12) then
489 : : if (abs(ig1min)<abs(ig1max)) ig1=abs(ig1max)
490 : : if (abs(ig1min)>abs(ig1max)) ig1=n1-abs(ig1min)
491 : : end if
492 : : if (abs(abs(qpt_(2))-half)<tol12) then
493 : : if (abs(ig2min)<abs(ig2max)) ig2=abs(ig2max)
494 : : if (abs(ig2min)>abs(ig2max)) ig2=n2-abs(ig2min)
495 : : end if
496 : : if (abs(abs(qpt_(3))-half)<tol12) then
497 : : if (abs(ig3min)<abs(ig3max)) ig3=abs(ig3max)
498 : : if (abs(ig3min)>abs(ig3max)) ig3=n3-abs(ig3min)
499 : : end if
500 : : call zerosym(work1,2,n1,n2,n3,ig1=ig1,ig2=ig2,ig3=ig3,&
501 : : comm_fft=mpi_enreg%comm_fft,distribfft=mpi_enreg%distribfft)
502 : : end if
503 : : end if
504 : :
505 : : ! Fourier Transform
506 : 0 : call fourdp(cplex,work1,vunitpawspherepot,1,mpi_enreg,nfft,1,ngfft,0)
507 : :
508 [ # # ]: 0 : ABI_FREE(work1)
509 : :
510 [ # # ]: 0 : end subroutine mkunitpawspherepot
511 : : !!***
512 : :
513 : :
514 : : !!****f* m_spacepar/hartre
515 : : !! NAME
516 : : !! hartre
517 : : !!
518 : : !! FUNCTION
519 : : !! Given rho(G), compute Hartree potential (=FFT of rho(G)/pi/(G+q)**2)
520 : : !! When cplex=1, assume q=(0 0 0), and vhartr will be REAL
521 : : !! When cplex=2, q must be taken into account, and vhartr will be COMPLEX
522 : : !!
523 : : !! NOTES
524 : : !! *Modified code to avoid if statements inside loops to skip G=0.
525 : : !! Replaced if statement on G^2>gsqcut to skip G s outside where
526 : : !! rho(G) should be 0. Effect is negligible but gsqcut should be
527 : : !! used to be strictly consistent with usage elsewhere in code.
528 : : !! *The speed-up is provided by doing a few precomputations outside
529 : : !! the inner loop. One variable size array is needed for this (gq).
530 : : !!
531 : : !! INPUTS
532 : : !! cplex= if 1, vhartr is REAL, if 2, vhartr is COMPLEX
533 : : !! gsqcut=cutoff value on G**2 for sphere inside fft box.
534 : : !! (gsqcut=(boxcut**2)*ecut/(2.d0*(Pi**2))
535 : : !! icutcoul= type of Coulomb cutoff to apply
536 : : !! izero=if 1, unbalanced components of Vhartree(g) are set to zero
537 : : !! mpi_enreg=information about MPI parallelization
538 : : !! nfft=(effective) number of FFT grid points (for this processor)
539 : : !! ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft
540 : : !! [qpt(3)=reduced coordinates for a wavevector to be combined with the G vectors (needed if cplex==2).]
541 : : !! rhog(2,nfft)=electron density in G space
542 : : !! rprimd(3,3)=dimensional primitive translations in real space (bohr)
543 : : !! vcutgeo(3)= array to describe the geometry of the Coulomb cutoff
544 : : !!
545 : : !! OUTPUT
546 : : !! vhartr(cplex*nfft)=Hartree potential in real space, either REAL or COMPLEX
547 : : !!
548 : : !! SOURCE
549 : :
550 : 84705 : subroutine hartre(cplex,gsqcut,icutcoul,izero,mpi_enreg,nfft,ngfft,nkpt,&
551 : 84705 : rcut,rhog,rprimd,vcutgeo,vhartr,&
552 : : qpt) ! Optional arguments
553 : :
554 : : !Arguments ------------------------------------
555 : : !scalars
556 : : integer,intent(in) :: cplex,icutcoul,izero,nfft,nkpt
557 : : real(dp),intent(in) :: gsqcut,rcut
558 : : type(MPI_type),intent(in) :: mpi_enreg
559 : : !arrays
560 : : integer,intent(in) :: ngfft(18)
561 : : real(dp),intent(in) :: rprimd(3,3),rhog(2,nfft),vcutgeo(3)
562 : : real(dp),intent(in),optional :: qpt(3)
563 : : real(dp),intent(out) :: vhartr(cplex*nfft)
564 : :
565 : : !Local variables-------------------------------
566 : : !scalars
567 : : integer,parameter :: im=2,re=1
568 : : integer :: i1,i2,i23,i2_local,i3,id1,id2,id3
569 : : integer :: ig,ig1min,ig1,ig1max,ig2,ig2min,ig2max,ig3,ig3min,ig3max
570 : : integer :: ii,ii1,ing,n1,n2,n3,qeq0,qeq05,me_fft,nproc_fft
571 : : integer :: nog0
572 : : real(dp),parameter :: tolfix=1.000000001e0_dp
573 : : real(dp) :: cutoff,den,gqg2p3,gqgm12,gqgm13,gqgm23,gs,gs2,gs3,ucvol
574 : : character(len=500) :: message
575 : : !arrays
576 : : integer :: id(3)
577 : 84705 : integer, contiguous, pointer :: fftn2_distrib(:),ffti2_local(:)
578 : 84705 : integer, contiguous, pointer :: fftn3_distrib(:),ffti3_local(:)
579 : : real(dp) :: gmet(3,3),gprimd(3,3),qpt_(3),rmet(3,3),tsec(2)
580 : 84705 : real(dp),allocatable :: gcutoff(:)
581 : 84705 : real(dp),allocatable :: gq(:,:),work1(:,:)
582 : : ! *************************************************************************
583 : :
584 : : ! Keep track of total time spent in hartre
585 : 84705 : call timab(10,1,tsec)
586 : :
587 : : ! Check that cplex has an allowed value
588 [ - + ]: 84705 : if (cplex/=1 .and. cplex/=2) then
589 : : write(message, '(a,i0,a,a)' )&
590 : 0 : 'From the calling routine, cplex=',cplex,ch10,'but the only value allowed are 1 and 2.'
591 : 0 : ABI_BUG(message)
592 : : end if
593 : :
594 : 84705 : call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
595 : :
596 : 84705 : n1=ngfft(1); n2=ngfft(2); n3=ngfft(3)
597 : 84705 : nproc_fft = mpi_enreg%nproc_fft; me_fft = mpi_enreg%me_fft
598 : :
599 : : ! Get the distrib associated with this fft_grid
600 : 84705 : call ptabs_fourdp(mpi_enreg,n2,n3,fftn2_distrib,ffti2_local,fftn3_distrib,ffti3_local)
601 : :
602 : : ! Initialize a few quantities
603 : 84705 : cutoff=gsqcut*tolfix
604 [ + + ]: 84705 : if(present(qpt))then
605 : 21140 : qpt_=qpt
606 : : else
607 : 63565 : qpt_=zero
608 : : end if
609 : 84705 : qeq0=0
610 [ + + ]: 84705 : if(qpt_(1)**2+qpt_(2)**2+qpt_(3)**2<1.d-15) qeq0=1
611 : 6064 : qeq05=0
612 : : if (qeq0==0) then
613 [ + + + + : 6064 : if (abs(abs(qpt_(1))-half)<tol12.or.abs(abs(qpt_(2))-half)<tol12.or.abs(abs(qpt_(3))-half)<tol12) qeq05=1
+ + ]
614 : : end if
615 : :
616 : : ! If cplex=1 then qpt_ should be 0 0 0
617 [ + + - + ]: 84705 : if (cplex==1.and. qeq0/=1) then
618 : : write(message,'(a,3e12.4,a,a)')&
619 : 0 : 'cplex=1 but qpt=',qpt_,ch10,&
620 : 0 : 'qpt should be 0 0 0.'
621 : 0 : ABI_BUG(message)
622 : : end if
623 : :
624 : : ! If FFT parallelism then qpt should not be 1/2
625 [ - + ]: 84705 : if (nproc_fft>1.and.qeq05==1) then
626 : : write(message, '(a,3e12.4,a,a)' )&
627 : 0 : 'FFT parallelism selected but qpt',qpt_,ch10,&
628 : 0 : 'qpt(i) should not be 1/2...'
629 : 0 : ABI_ERROR(message)
630 : : end if
631 : :
632 : : !PCM cut-off is implemented outside termcutoff
633 [ + + + + ]: 84705 : nog0=0; if (qeq0==1 .or. icutcoul==55) nog0=1
634 : :
635 : : !Initialize Gcut-off array from m_gtermcutoff
636 : 84705 : call termcutoff(gcutoff,gsqcut,icutcoul,ngfft,nkpt,rcut,rprimd,vcutgeo,qpt=qpt_)
637 : :
638 : : ! In order to speed the routine, precompute the components of g+q
639 : : ! Also check if the booked space was large enough...
640 [ - + + - : 254115 : ABI_MALLOC(gq,(3,max(n1,n2,n3)))
+ - + - -
+ + - ]
641 [ + + ]: 338820 : do ii=1,3
642 : 254115 : id(ii)=ngfft(ii)/2+2
643 [ + + ]: 6105429 : do ing=1,ngfft(ii)
644 : 5766609 : ig=ing-(ing/id(ii))*ngfft(ii)-1
645 : 6020724 : gq(ii,ing)=ig+qpt_(ii)
646 : : end do
647 : : end do
648 : 84705 : ig1max=-1;ig2max=-1;ig3max=-1
649 : 84705 : ig1min=n1;ig2min=n2;ig3min=n3
650 : :
651 [ - + + - : 254115 : ABI_MALLOC(work1,(2,nfft))
+ - + - -
+ - + ]
652 : 84705 : id1=n1/2+2;id2=n2/2+2;id3=n3/2+2
653 : :
654 : : ! Triple loop on each dimension
655 [ + + ]: 2061476 : do i3=1,n3
656 : 1976771 : ig3=i3-(i3/id3)*n3-1
657 : : ! Precompute some products that do not depend on i2 and i1
658 : 1976771 : gs3=gq(3,i3)*gq(3,i3)*gmet(3,3)
659 : 1976771 : gqgm23=gq(3,i3)*gmet(2,3)*2
660 : 1976771 : gqgm13=gq(3,i3)*gmet(1,3)*2
661 : :
662 [ + + ]: 55754389 : do i2=1,n2
663 : 53692913 : ig2=i2-(i2/id2)*n2-1
664 [ + + ]: 55669684 : if (fftn2_distrib(i2) == me_fft) then
665 : 51359633 : gs2=gs3+ gq(2,i2)*(gq(2,i2)*gmet(2,2)+gqgm23)
666 : 51359633 : gqgm12=gq(2,i2)*gmet(1,2)*2
667 : 51359633 : gqg2p3=gqgm13+gqgm12
668 : :
669 : 51359633 : i2_local = ffti2_local(i2)
670 : 51359633 : i23=n1*(i2_local-1 +(n2/nproc_fft)*(i3-1))
671 : : ! Do the test that eliminates the Gamma point outside of the inner loop
672 : 51359633 : ii1=1
673 : : ! if(i23==0 .and. qeq0==1 .and. ig2==0 .and. ig3==0)then
674 [ + + + + ]: 51359633 : if(i23==0 .and. nog0==1 .and. ig2==0 .and. ig3==0)then
675 : 76845 : ii1=2
676 : 76845 : work1(re,1+i23)=zero
677 : 76845 : work1(im,1+i23)=zero
678 : : end if
679 : :
680 : : ! Final inner loop on the first dimension (note the lower limit)
681 [ + + ]: 2105889925 : do i1=ii1,n1
682 : 2054530292 : gs=gs2+ gq(1,i1)*(gq(1,i1)*gmet(1,1)+gqg2p3)
683 : 2054530292 : ii=i1+i23
684 : :
685 [ + + ]: 2105889925 : if(gs<=cutoff)then
686 : : ! Identify min/max indexes (to cancel unbalanced contributions later)
687 : : ! Count (q+g)-vectors with similar norm
688 [ + + - + ]: 807379246 : if ((qeq05==1).and.(izero==1)) then
689 : 0 : ig1=i1-(i1/id1)*n1-1
690 : 0 : ig1max=max(ig1max,ig1); ig1min=min(ig1min,ig1)
691 : 0 : ig2max=max(ig2max,ig2); ig2min=min(ig2min,ig2)
692 : 0 : ig3max=max(ig3max,ig3); ig3min=min(ig3min,ig3)
693 : : end if
694 : :
695 : 807379246 : den=piinv/gs*gcutoff(ii)
696 : 807379246 : work1(re,ii)=rhog(re,ii)*den
697 : 807379246 : work1(im,ii)=rhog(im,ii)*den
698 : : else
699 : : ! gs>cutoff
700 : 1247151046 : work1(re,ii)=zero
701 : 1247151046 : work1(im,ii)=zero
702 : : end if
703 : :
704 : : end do ! End loop on i1
705 : : end if
706 : : end do ! End loop on i2
707 : : end do ! End loop on i3
708 : :
709 [ - + ]: 84705 : ABI_FREE(gq)
710 : :
711 [ + + ]: 84705 : if (izero==1) then
712 : : ! Set contribution of unbalanced components to zero
713 : :
714 [ + - ]: 23896 : if (qeq0==1) then !q=0
715 : 23896 : call zerosym(work1,2,n1,n2,n3,comm_fft=mpi_enreg%comm_fft,distribfft=mpi_enreg%distribfft)
716 : :
717 [ # # ]: 0 : else if (qeq05==1) then
718 : : !q=1/2; this doesn't work in parallel
719 [ # # ]: 0 : ig1=-1;if (mod(n1,2)==0) ig1=1+n1/2
720 [ # # ]: 0 : ig2=-1;if (mod(n2,2)==0) ig2=1+n2/2
721 [ # # ]: 0 : ig3=-1;if (mod(n3,2)==0) ig3=1+n3/2
722 [ # # ]: 0 : if (abs(abs(qpt_(1))-half)<tol12) then
723 [ # # ]: 0 : if (abs(ig1min)<abs(ig1max)) ig1=abs(ig1max)
724 [ # # ]: 0 : if (abs(ig1min)>abs(ig1max)) ig1=n1-abs(ig1min)
725 : : end if
726 [ # # ]: 0 : if (abs(abs(qpt_(2))-half)<tol12) then
727 [ # # ]: 0 : if (abs(ig2min)<abs(ig2max)) ig2=abs(ig2max)
728 [ # # ]: 0 : if (abs(ig2min)>abs(ig2max)) ig2=n2-abs(ig2min)
729 : : end if
730 [ # # ]: 0 : if (abs(abs(qpt_(3))-half)<tol12) then
731 [ # # ]: 0 : if (abs(ig3min)<abs(ig3max)) ig3=abs(ig3max)
732 [ # # ]: 0 : if (abs(ig3min)>abs(ig3max)) ig3=n3-abs(ig3min)
733 : : end if
734 : : call zerosym(work1,2,n1,n2,n3,ig1=ig1,ig2=ig2,ig3=ig3,&
735 : 0 : comm_fft=mpi_enreg%comm_fft,distribfft=mpi_enreg%distribfft)
736 : : end if
737 : : end if
738 : :
739 : : ! Fourier Transform Vhartree. Vh in reciprocal space was stored in work1
740 : 84705 : call fourdp(cplex,work1,vhartr,1,mpi_enreg,nfft,1,ngfft,0)
741 : :
742 [ - + ]: 84705 : ABI_FREE(gcutoff)
743 [ - + ]: 84705 : ABI_FREE(work1)
744 : :
745 : 84705 : call timab(10,2,tsec)
746 : :
747 : 84705 : end subroutine hartre
748 : : !!***
749 : :
750 : : !!****f* m_spacepar/meanvalue_g
751 : : !! NAME
752 : : !! meanvalue_g
753 : : !!
754 : : !! FUNCTION
755 : : !! Compute the mean value of one wavefunction, in reciprocal space,
756 : : !! for an operator that is real, diagonal in G-space: <wf|op|wf>
757 : : !! For the time being, only spin-independent operators are treated.
758 : : !!
759 : : !! INPUTS
760 : : !! diag(npw)=diagonal operator (real, spin-independent!)
761 : : !! filter= if 1, need to filter on the value of diag, that must be less than huge(zero)*1.d-11
762 : : !! otherwise, should be 0
763 : : !! istwf_k=storage mode of the vectors
764 : : !! npw=number of planewaves of the vector
765 : : !! nspinor=number of spinor components
766 : : !! vect(2,npw*nspinor)=vector
767 : : !! vect1(2,npw*nspinor)=vector1 (=vector in most of the cases)
768 : : !! use_ndo = says if vect=/vect1
769 : : !!
770 : : !! OUTPUT
771 : : !! ar=mean value
772 : : !!
773 : : !! SOURCE
774 : :
775 : 5390029 : subroutine meanvalue_g(ar,diag,filter,istwf_k,mpi_enreg,npw,nspinor,vect,vect1,use_ndo,ar_im, &
776 : : gpu_thread_limit) ! optional
777 : :
778 : : !Arguments ------------------------------------
779 : : !scalars
780 : : integer,intent(in) :: filter,istwf_k,npw,nspinor,use_ndo
781 : : integer,intent(in),optional :: gpu_thread_limit
782 : : real(dp),intent(out) :: ar
783 : : real(dp),intent(out),optional :: ar_im
784 : : type(MPI_type),intent(in) :: mpi_enreg
785 : : !arrays
786 : : real(dp),intent(in) :: diag(npw),vect(2,npw*nspinor)
787 : : real(dp),intent(in) :: vect1(2,npw*nspinor)
788 : :
789 : : !Local variables-------------------------------
790 : : !scalars
791 : : integer :: i1,ierr,ipw,jpw,me_g0,nthreads_bak,l_gpu_thread_limit
792 : : character(len=500) :: message
793 : : real(dp), parameter :: hugevalue = huge(zero)*1.d-11
794 : : ! *************************************************************************
795 : :
796 : : DBG_CHECK(ANY(filter==(/0,1/)),"Wrong filter")
797 : : DBG_CHECK(ANY(nspinor==(/1,2/)),"Wrong nspinor")
798 : : DBG_CHECK(ANY(istwf_k==(/(ipw,ipw=1,9)/)),"Wrong istwf_k")
799 : :
800 [ + + - + ]: 5390029 : if(nspinor==2 .and. istwf_k/=1)then
801 : : write(message,'(3a,i0,a,i0)')&
802 : 0 : 'When istwf_k/=1, nspinor must be 1,',ch10, 'however, nspinor=',nspinor,', and istwf_k=',istwf_k
803 : 0 : ABI_BUG(message)
804 : : end if
805 : :
806 [ + + - + : 5390029 : if(use_ndo==1 .and. (istwf_k==2 .and.me_g0==1)) then
- - ]
807 : 0 : ABI_BUG('use_ndo==1, not tested, use istwfk=1')
808 : : end if
809 : :
810 : 5390029 : me_g0 = mpi_enreg%me_g0
811 : :
812 : 5390029 : ar=zero
813 [ + + ]: 5390029 : if(present(ar_im)) ar_im=zero
814 : :
815 [ + + ]: 5390029 : l_gpu_thread_limit=0; if(present(gpu_thread_limit)) l_gpu_thread_limit=gpu_thread_limit
816 : :
817 [ - + ]: 4209085 : if(l_gpu_thread_limit /= 0) then
818 : 0 : nthreads_bak=xomp_get_max_threads()
819 : 0 : call xomp_set_num_threads(min(l_gpu_thread_limit,nthreads_bak))
820 : : end if
821 : :
822 : : !Normal storage mode
823 [ + + ]: 5390029 : if(istwf_k==1)then
824 : :
825 : : ! No filter
826 [ + + ]: 5219677 : if(filter==0)then
827 : : !$OMP PARALLEL DO REDUCTION(+:ar)
828 [ + + ]: 1666261859 : do ipw=1,npw
829 : 1666261859 : ar=ar+diag(ipw)*(vect(1,ipw)*vect1(1,ipw)+vect(2,ipw)*vect1(2,ipw))
830 : : end do
831 [ + + ]: 5144541 : if(nspinor==2)then
832 : : !$OMP PARALLEL DO REDUCTION(+:ar) PRIVATE(jpw)
833 [ + + ]: 716684132 : do ipw=1+npw,2*npw
834 : 716150434 : jpw=ipw-npw
835 : 716684132 : ar=ar+diag(jpw)*(vect(1,ipw)*vect1(1,ipw)+vect(2,ipw)*vect1(2,ipw))
836 : : end do
837 : : end if
838 [ + + ]: 5144541 : if(use_ndo==1)then
839 : : !$OMP PARALLEL DO REDUCTION(+:ar_im)
840 [ + + ]: 825657098 : do ipw=1,npw
841 : 825657098 : ar_im=ar_im+diag(ipw)*(vect1(1,ipw)*vect(2,ipw)-vect1(2,ipw)*vect(1,ipw))
842 : : end do
843 [ + + ]: 756112 : if(nspinor == 2) then
844 : : !$OMP PARALLEL DO REDUCTION(+:ar_im) PRIVATE(jpw)
845 [ + + ]: 634851412 : do ipw=1+npw,2*npw
846 : 634601868 : jpw=ipw-npw
847 : 634851412 : ar_im=ar_im+diag(jpw)*(vect1(1,ipw)*vect(2,ipw)-vect1(2,ipw)*vect(1,ipw))
848 : : end do
849 : : end if
850 : : end if
851 : :
852 : : else
853 : : ! will filter
854 : :
855 : : !$OMP PARALLEL DO REDUCTION(+:ar)
856 [ + + ]: 26463664 : do ipw=1,npw
857 [ + + ]: 26463664 : if(diag(ipw)<hugevalue)then
858 : 25769392 : ar=ar+diag(ipw)*(vect(1,ipw)*vect1(1,ipw)+vect(2,ipw)*vect1(2,ipw))
859 : : end if
860 : : end do
861 [ - + ]: 75136 : if(nspinor==2)then
862 : : !$OMP PARALLEL DO REDUCTION(+:ar) PRIVATE(jpw)
863 [ # # ]: 0 : do ipw=1+npw,2*npw
864 : 0 : jpw=ipw-npw
865 [ # # ]: 0 : if(diag(jpw)<hugevalue)then
866 : 0 : ar=ar+diag(jpw)*(vect(1,ipw)*vect1(1,ipw)+vect(2,ipw)*vect1(2,ipw))
867 : : end if
868 : : end do
869 : : end if
870 [ - + ]: 75136 : if(use_ndo==1)then
871 [ # # ]: 0 : if(.not.present(ar_im)) then
872 : 0 : ABI_BUG("use_ndo true and ar_im not present")
873 : : end if
874 : : !$OMP PARALLEL DO REDUCTION(+:ar_im)
875 [ # # ]: 0 : do ipw=1,npw
876 [ # # ]: 0 : if(diag(ipw)<hugevalue)then
877 : 0 : ar_im=ar_im+diag(ipw)*(vect1(1,ipw)*vect(2,ipw)-vect1(2,ipw)*vect(1,ipw))
878 : : end if
879 : : end do
880 [ # # ]: 0 : if(nspinor == 2) then
881 : : !$OMP PARALLEL DO REDUCTION(+:ar_im) PRIVATE(jpw)
882 [ # # ]: 0 : do ipw=1+npw,2*npw
883 : 0 : jpw=ipw-npw
884 [ # # ]: 0 : if(diag(jpw)<hugevalue)then
885 : 0 : ar_im=ar_im+diag(jpw)*(vect1(1,ipw)*vect(2,ipw)-vect1(2,ipw)*vect(1,ipw))
886 : : end if
887 : : end do
888 : : end if
889 : : end if
890 : :
891 : : end if ! filter==0
892 : :
893 [ + - ]: 170352 : else if(istwf_k>=2)then
894 : :
895 [ + - ]: 170352 : if(filter==0)then
896 : 170352 : i1=1
897 [ + + + + ]: 170352 : if(istwf_k==2 .and. me_g0==1)then ! MPIWF need to know which proc has G=0
898 : 73317 : ar=half*diag(1)*vect(1,1)*vect1(1,1) ; i1=2
899 : : end if
900 : :
901 : : !$OMP PARALLEL DO REDUCTION(+:ar)
902 [ + + ]: 53195775 : do ipw=i1,npw
903 : 53195775 : ar=ar+diag(ipw)*(vect(1,ipw)*vect1(1,ipw)+vect(2,ipw)*vect1(2,ipw))
904 : : end do
905 : :
906 : :
907 : : else ! filter/=0
908 : 0 : i1=1
909 [ # # # # ]: 0 : if(istwf_k==2 .and. me_g0==1)then
910 [ # # ]: 0 : if(diag(1)<hugevalue)then
911 : 0 : ar=half*diag(1)*vect(1,1)*vect1(1,1) ; i1=2
912 : : end if
913 : : end if
914 : :
915 : : !$OMP PARALLEL DO REDUCTION(+:ar)
916 [ # # ]: 0 : do ipw=i1,npw
917 [ # # ]: 0 : if(diag(ipw)<hugevalue)then
918 : 0 : ar=ar+diag(ipw)*(vect(1,ipw)*vect1(1,ipw)+vect(2,ipw)*vect1(2,ipw))
919 : : end if
920 : : end do
921 : : end if ! filter==0
922 : :
923 : 170352 : ar=two*ar
924 : :
925 : : end if ! istwf_k
926 : :
927 : : !MPIWF need to make reduction on ar and ai .
928 [ + + ]: 5390029 : if(mpi_enreg%paral_kgb==1)then
929 : 427100 : call xmpi_sum(ar,mpi_enreg%comm_bandspinorfft ,ierr)
930 [ + + ]: 427100 : if (present(ar_im)) call xmpi_sum(ar_im,mpi_enreg%comm_bandspinorfft,ierr)
931 : : end if
932 : :
933 [ - + ]: 5390029 : if (l_gpu_thread_limit /= 0) call xomp_set_num_threads(nthreads_bak)
934 : :
935 : 5390029 : end subroutine meanvalue_g
936 : : !!***
937 : :
938 : : !!****f* m_spacepar/meanvalue_g_batch
939 : : !! NAME
940 : : !! meanvalue_g_batch
941 : : !!
942 : : !! FUNCTION
943 : : !! Batched version of meanvalue_g: computes ndat mean values <psi_i|op|psi_i>
944 : : !! for wavefunctions packed contiguously in a single array, where op is real
945 : : !! and diagonal in G-space.
946 : : !!
947 : : !! Note: nspinor must be 1 when istwf_k/=1 (same constraint as meanvalue_g).
948 : : !! Other combinations fall back to scalar meanvalue_g calls.
949 : : !!
950 : : !! INPUTS
951 : : !! diag(npw)=diagonal operator (real, spin-independent)
952 : : !! filter= if 1, filter on diag < huge*1.d-11; otherwise 0
953 : : !! istwf_k=storage mode of the vectors
954 : : !! npw=number of planewaves per wavefunction
955 : : !! nspinor=number of spinor components
956 : : !! ndat=number of wavefunctions (batch size)
957 : : !! vect(2,npw*nspinor*ndat)=packed input wavefunctions; band idat occupies
958 : : !! columns 1+(idat-1)*npw*nspinor : idat*npw*nspinor
959 : : !! vect1(2,npw*nspinor*ndat)=second set of wavefunctions (equals vect when use_ndo==0)
960 : : !! use_ndo=1 if vect /= vect1 (non-diagonal operator); 0 otherwise
961 : : !! gpu_option= (optional) GPU acceleration flag (ABI_GPU_OPENMP, etc.)
962 : : !!
963 : : !! OUTPUT
964 : : !! ar(ndat)=mean values, one per wavefunction
965 : : !!
966 : : !! SOURCE
967 : :
968 : 2129216 : subroutine meanvalue_g_batch(ar, diag, filter, istwf_k, mpi_enreg, npw, nspinor, ndat, &
969 : 2129216 : vect, vect1, use_ndo,&
970 : : gpu_option, gpu_thread_limit) ! optional
971 : :
972 : : !Arguments ------------------------------------
973 : : !scalars
974 : : integer, intent(in) :: filter, istwf_k, npw, nspinor, ndat, use_ndo
975 : : integer, intent(in), optional :: gpu_option, gpu_thread_limit
976 : : type(MPI_type), intent(in) :: mpi_enreg
977 : : !arrays
978 : : real(dp), intent(out) :: ar(ndat)
979 : : real(dp), target, intent(in) :: diag(npw)
980 : : real(dp), target, intent(in) :: vect(2, npw*nspinor*ndat)
981 : : real(dp), target, intent(in) :: vect1(2, npw*nspinor*ndat)
982 : :
983 : : !Local variables-------------------------------
984 : : !scalars
985 : : integer :: idat, ipw, jpw, ierr, l_gpu_option, l_gpu_thread_limit, nthreads_bak, i1, me_g0
986 : : real(dp) :: local_ar
987 : : character(len=500) :: message
988 : : ! *************************************************************************
989 : :
990 [ + + - + ]: 2129216 : if(nspinor==2 .and. istwf_k/=1)then
991 : : write(message,'(a,a,a,i6,a,i6)')&
992 : 0 : 'When istwf_k/=1, nspinor must be 1,',ch10,&
993 : 0 : 'however, nspinor=',nspinor,', and istwf_k=',istwf_k
994 : 0 : ABI_BUG(message)
995 : : end if
996 : :
997 [ + - ]: 2129216 : l_gpu_option = ABI_GPU_DISABLED; if (present(gpu_option)) l_gpu_option = gpu_option
998 [ - + ]: 2129216 : l_gpu_thread_limit=0; if(present(gpu_thread_limit)) l_gpu_thread_limit=gpu_thread_limit
999 [ - + ]: 2129216 : if(l_gpu_option==ABI_GPU_OPENMP) l_gpu_thread_limit=0
1000 : :
1001 [ + + ]: 4734780 : ar(:) = zero
1002 : :
1003 [ + + + - : 2129216 : if (istwf_k == 1 .and. filter == 0 .and. use_ndo == 0) then
+ - ]
1004 : :
1005 : : #ifdef HAVE_OPENMP_OFFLOAD
1006 : : !$OMP TARGET TEAMS DISTRIBUTE MAP(to:diag,vect) MAP(tofrom:ar) &
1007 : : !$OMP& IF(l_gpu_option==ABI_GPU_OPENMP)
1008 : : #endif
1009 [ + + ]: 4123940 : do idat = 1, ndat
1010 : : local_ar = zero
1011 : : !$OMP PARALLEL DO REDUCTION(+:local_ar) PRIVATE(jpw)
1012 [ + + ]: 657889592 : do ipw = 1, npw*nspinor
1013 : 655686129 : jpw = mod(ipw-1, npw) + 1
1014 : : local_ar = local_ar + diag(jpw) * ( vect(1, ipw+(idat-1)*npw*nspinor)**2 &
1015 : 657889592 : +vect(2, ipw+(idat-1)*npw*nspinor)**2)
1016 : : end do
1017 : 4123940 : ar(idat) = local_ar
1018 : : end do
1019 : :
1020 [ + + ]: 1920477 : if (mpi_enreg%paral_kgb == 1) then
1021 : 112328 : call xmpi_sum(ar, ndat, mpi_enreg%comm_bandspinorfft, ierr)
1022 : : end if
1023 : :
1024 [ - + - + : 208739 : else if (istwf_k >= 2 .and. filter == 0 .and. use_ndo == 0) then
- + ]
1025 : :
1026 : 208739 : me_g0 = mpi_enreg%me_g0
1027 : 208739 : i1 = 1
1028 [ + + + + ]: 208739 : if (istwf_k == 2 .and. me_g0 == 1) i1 = 2
1029 : :
1030 : : #ifdef HAVE_OPENMP_OFFLOAD
1031 : : !$OMP TARGET TEAMS DISTRIBUTE MAP(to:diag,vect) MAP(tofrom:ar) &
1032 : : !$OMP& IF(l_gpu_option==ABI_GPU_OPENMP)
1033 : : #endif
1034 [ + + ]: 610840 : do idat = 1, ndat
1035 : 402101 : local_ar = zero
1036 [ + + ]: 402101 : if (i1 == 2) then
1037 : 114825 : local_ar = half * diag(1) * vect(1, 1+(idat-1)*npw)**2
1038 : : end if
1039 : : !$OMP PARALLEL DO REDUCTION(+:local_ar)
1040 [ + + ]: 149034146 : do ipw = i1, npw
1041 : : local_ar = local_ar + diag(ipw) * (vect(1, ipw+(idat-1)*npw)**2 &
1042 : 149034146 : +vect(2, ipw+(idat-1)*npw)**2)
1043 : : end do
1044 : 610840 : ar(idat) = two * local_ar
1045 : : end do
1046 : :
1047 [ + + ]: 208739 : if (mpi_enreg%paral_kgb == 1) then
1048 : 21532 : call xmpi_sum(ar, ndat, mpi_enreg%comm_bandspinorfft, ierr)
1049 : : end if
1050 : :
1051 : : else
1052 : : ! Fallback: scalar loop for cases not yet GPU-ported (filter==1, use_ndo==1).
1053 : : ! Each call does its own MPI reduction.
1054 : : #ifdef HAVE_OPENMP_OFFLOAD
1055 : : if(l_gpu_option==ABI_GPU_OPENMP) then
1056 : : if(xomp_target_is_present(c_loc(diag))) then
1057 : : !$OMP TARGET UPDATE FROM(diag)
1058 : : end if
1059 : : if(xomp_target_is_present(c_loc(vect))) then
1060 : : !$OMP TARGET UPDATE FROM(vect)
1061 : : end if
1062 : : if(xomp_target_is_present(c_loc(vect1)) .and. .not. c_associated(c_loc(vect1), c_loc(vect))) then
1063 : : !$OMP TARGET UPDATE FROM(vect1)
1064 : : end if
1065 : : end if
1066 : : #endif
1067 [ # # ]: 0 : do idat = 1, ndat
1068 : : call meanvalue_g(ar(idat), diag, filter, istwf_k, mpi_enreg, npw, nspinor, &
1069 : : vect (:, 1+(idat-1)*npw*nspinor:idat*npw*nspinor), &
1070 : : vect1(:, 1+(idat-1)*npw*nspinor:idat*npw*nspinor), &
1071 : 0 : use_ndo)
1072 : : end do
1073 : :
1074 : : end if
1075 [ - + ]: 2129216 : if (l_gpu_thread_limit /= 0) call xomp_set_num_threads(nthreads_bak)
1076 : :
1077 : 2129216 : end subroutine meanvalue_g_batch
1078 : : !!***
1079 : :
1080 : : !!****f* m_spacepar/laplacian
1081 : : !! NAME
1082 : : !! laplacian
1083 : : !!
1084 : : !! FUNCTION
1085 : : !! compute the laplacian of a function defined in real space
1086 : : !! the code is written in the way of /3xc/xcden.F90
1087 : : !!
1088 : : !! INPUTS
1089 : : !! gprimd(3,3)=dimensional reciprocal space primitive translations
1090 : : !! mpi_enreg=information about MPI parallelization
1091 : : !! nfft=number of points of the fft grid
1092 : : !! nfunc=number of functions on the grid for which the laplacian is to be calculated
1093 : : !! ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft
1094 : : !! (optional) rdfuncr(nfft,nfunc)=real(dp) discretized functions in real space
1095 : : !! rdfuncg_in TO BE DESCRIBED SB 090901
1096 : : !! laplacerdfuncg_in TO BE DESCRIBED SB 090901
1097 : : !! (optional) g2cart_in(nfft) = G**2 on the grid
1098 : : !!
1099 : : !! OUTPUT
1100 : : !! (optional) laplacerdfuncr = laplacian in real space of the functions in rdfuncr
1101 : : !! (optional) rdfuncg = real(dp) discretized functions in fourier space
1102 : : !! (optional) laplacerdfuncg = real(dp) discretized laplacian of the functions in fourier space
1103 : : !! (optional) g2cart_out(nfft) = G**2 on the grid
1104 : : !! rdfuncg_out TO BE DESCRIBED SB 090901
1105 : : !! laplacerdfuncg_out TO BE DESCRIBED SB 090901
1106 : : !!
1107 : : !! SOURCE
1108 : :
1109 : 8148 : subroutine laplacian(gprimd,mpi_enreg,nfft,nfunc,ngfft,rdfuncr,&
1110 : 8148 : laplacerdfuncr,rdfuncg_out,laplacerdfuncg_out,g2cart_out,rdfuncg_in,g2cart_in)
1111 : :
1112 : : !Arguments ------------------------------------
1113 : : !scalars
1114 : : integer,intent(in) :: nfft,nfunc
1115 : : type(MPI_type),intent(in) :: mpi_enreg
1116 : : !arrays
1117 : : integer,intent(in) :: ngfft(18)
1118 : : real(dp),intent(in) :: gprimd(3,3)
1119 : : real(dp),intent(inout),optional :: laplacerdfuncr(nfft,nfunc)
1120 : : real(dp),intent(inout),optional,target :: rdfuncr(nfft,nfunc)
1121 : : real(dp),intent(in),optional,target :: g2cart_in(nfft) !vz_i
1122 : : real(dp),intent(out),optional,target :: g2cart_out(nfft) !vz_i
1123 : : real(dp),intent(out),optional,target :: laplacerdfuncg_out(2,nfft,nfunc)
1124 : : real(dp),intent(in),optional,target :: rdfuncg_in(2,nfft,nfunc) !vz_i
1125 : : real(dp),intent(out),optional,target :: rdfuncg_out(2,nfft,nfunc)
1126 : :
1127 : : !Local variables-------------------------------
1128 : : !scalars
1129 : : integer :: count,i1,i2,i3,id1,id2,id3,ifft,ifunc,ig1,ig2,ig3,ii1,n1,n2
1130 : : integer :: n3
1131 : : real(dp) :: b11,b12,b13,b21,b22,b23,b31,b32,b33
1132 : : !arrays
1133 : 8148 : integer, contiguous, pointer :: fftn2_distrib(:),ffti2_local(:)
1134 : 8148 : integer, contiguous, pointer :: fftn3_distrib(:),ffti3_local(:)
1135 : 8148 : real(dp),contiguous, pointer :: g2cart(:),laplacerdfuncg(:,:,:),rdfuncg(:,:,:)
1136 : : ! *************************************************************************
1137 : :
1138 : : !Keep local copy of fft dimensions
1139 : 8148 : n1=ngfft(1)
1140 : 8148 : n2=ngfft(2)
1141 : 8148 : n3=ngfft(3)
1142 : :
1143 : 8148 : if(present(laplacerdfuncg_out)) then
1144 : 0 : laplacerdfuncg => laplacerdfuncg_out
1145 : : else
1146 [ - + + - : 32592 : ABI_MALLOC(laplacerdfuncg,(2,nfft,nfunc))
- + + - +
- + - - +
- + ]
1147 : : end if
1148 : :
1149 : : ! Get the distrib associated with this fft_grid
1150 : 8148 : call ptabs_fourdp(mpi_enreg,n2,n3,fftn2_distrib,ffti2_local,fftn3_distrib,ffti3_local)
1151 : :
1152 : : !change the real density rdfuncr on real space on the real density
1153 : : !rdfuncg in reciprocal space
1154 [ + - ]: 8148 : if(.not.present(rdfuncg_in)) then
1155 [ - + ]: 8148 : if(present(rdfuncg_out)) then
1156 : 0 : rdfuncg => rdfuncg_out
1157 : : else
1158 [ - + + - : 32592 : ABI_MALLOC(rdfuncg,(2,nfft,nfunc))
- + + - +
- + - - +
- + ]
1159 : : end if
1160 [ + - ]: 8148 : if(present(rdfuncr)) then
1161 [ + + ]: 16296 : do ifunc=1,nfunc
1162 : 16296 : call fourdp(1,rdfuncg(:,:,ifunc),rdfuncr(:,ifunc),-1,mpi_enreg,nfft,1,ngfft,0)
1163 : : end do
1164 : : end if
1165 : : else
1166 : 0 : rdfuncg => rdfuncg_in
1167 : : end if
1168 : :
1169 : : !apply the laplacian on laplacerdfuncr
1170 : : !code from /3xc/xcden.F90
1171 : : !see also src/5common/hatre.F90 and src/5common/moddiel.F90
1172 : : !Keep local copy of fft dimensions
1173 : : !Initialize computation of G^2 in cartesian coordinates
1174 [ + + ]: 8148 : if(.not.present(g2cart_in)) then
1175 [ + + ]: 2239 : if(present(g2cart_out)) then
1176 : 7 : g2cart => g2cart_out
1177 : : else
1178 [ - + + - : 6696 : ABI_MALLOC(g2cart,(nfft))
+ - + - -
+ - + ]
1179 : : end if
1180 : 2239 : id1=int(n1/2)+2
1181 : 2239 : id2=int(n2/2)+2
1182 : 2239 : id3=int(n3/2)+2
1183 : 2239 : count=0
1184 [ + + ]: 226139 : do i3=1,n3
1185 : 223900 : ifft=(i3-1)*n1*(n2/mpi_enreg%nproc_fft)
1186 : 223900 : ig3=i3-int(i3/id3)*n3-1
1187 [ + + ]: 2465139 : do i2=1,n2
1188 [ + - ]: 2462900 : if (fftn2_distrib(i2)==mpi_enreg%me_fft) then
1189 : 2239000 : ig2=i2-int(i2/id2)*n2-1
1190 : :
1191 : 2239000 : ii1=1
1192 [ + + ]: 24629000 : do i1=ii1,n1
1193 : 22390000 : ig1=i1-int(i1/id1)*n1-1
1194 : 22390000 : ifft=ifft+1
1195 : :
1196 : 22390000 : b11=gprimd(1,1)*real(ig1,dp)
1197 : 22390000 : b21=gprimd(2,1)*real(ig1,dp)
1198 : 22390000 : b31=gprimd(3,1)*real(ig1,dp)
1199 : 22390000 : b12=gprimd(1,2)*real(ig2,dp)
1200 : 22390000 : b22=gprimd(2,2)*real(ig2,dp)
1201 : 22390000 : b32=gprimd(3,2)*real(ig2,dp)
1202 : 22390000 : b13=gprimd(1,3)*real(ig3,dp)
1203 : 22390000 : b23=gprimd(2,3)*real(ig3,dp)
1204 : 22390000 : b33=gprimd(3,3)*real(ig3,dp)
1205 : :
1206 : : g2cart(ifft)=( &
1207 : : & (b11+b12+b13)**2&
1208 : : & +(b21+b22+b23)**2&
1209 : : & +(b31+b32+b33)**2&
1210 : 22390000 : & )
1211 [ + + ]: 47019000 : do ifunc=1,nfunc
1212 : : ! compute the laplacian in Fourier space that is * (i x 2pi x G)**2
1213 : 22390000 : laplacerdfuncg(1,ifft,ifunc) = -rdfuncg(1,ifft,ifunc)*g2cart(ifft)*two_pi*two_pi
1214 : 44780000 : laplacerdfuncg(2,ifft,ifunc) = -rdfuncg(2,ifft,ifunc)*g2cart(ifft)*two_pi*two_pi
1215 : : end do
1216 : : end do
1217 : : end if
1218 : : end do
1219 : : end do
1220 [ + + ]: 2239 : if(.not.present(g2cart_out)) then
1221 [ - + ]: 2232 : ABI_FREE(g2cart)
1222 : : end if
1223 : : else
1224 : 5909 : g2cart => g2cart_in
1225 [ + + ]: 11818 : do ifunc=1,nfunc
1226 [ + + ]: 59101818 : do ifft=1,nfft
1227 : : ! compute the laplacian in Fourier space that is * (i x 2pi x G)**2
1228 : 59090000 : laplacerdfuncg(1,ifft,ifunc) = -rdfuncg(1,ifft,ifunc)*g2cart(ifft)*two_pi*two_pi
1229 : 59095909 : laplacerdfuncg(2,ifft,ifunc) = -rdfuncg(2,ifft,ifunc)*g2cart(ifft)*two_pi*two_pi
1230 : : end do
1231 : : end do
1232 : : end if
1233 : :
1234 : : !get the result back into real space
1235 [ + - ]: 8148 : if(present(laplacerdfuncr)) then
1236 [ + + ]: 16296 : do ifunc=1,nfunc
1237 : 16296 : call fourdp(1,laplacerdfuncg(:,:,ifunc),laplacerdfuncr(:,ifunc),1,mpi_enreg,nfft,1,ngfft,0)
1238 : : end do
1239 : : end if
1240 : :
1241 : : !deallocate pointers
1242 [ + - ]: 8148 : if((.not.present(rdfuncg_in)).and.(.not.present(rdfuncg_in))) then
1243 [ - + ]: 8148 : ABI_FREE(rdfuncg)
1244 : : end if
1245 [ + - ]: 8148 : if(.not.present(laplacerdfuncg_out)) then
1246 [ - + ]: 8148 : ABI_FREE(laplacerdfuncg)
1247 : : end if
1248 : :
1249 [ - + - + : 24444 : end subroutine laplacian
- + + - +
- - + ]
1250 : : !!***
1251 : :
1252 : : !!****f* m_spacepar/redgr
1253 : : !! NAME
1254 : : !! redgr
1255 : : !!
1256 : : !! FUNCTION
1257 : : !! Compute reduced gradients of a real function on the usual unshifted
1258 : : !! fft grid. The gradient directions are the along the primitive
1259 : : !! reciprocal lattice vectors.
1260 : : !! The input function is intended to be a single spin component of
1261 : : !! the valence charge density, the valence + core charge densities
1262 : : !! or the first-order core charge density for use in frozen wf
1263 : : !! elastic tensor calculations within the GGA.
1264 : : !!
1265 : : !! NOTES
1266 : : !! Closely linked to xcden, but limited to Q=0, real charge densities,
1267 : : !! and unshifted grids.
1268 : : !!
1269 : : !! INPUTS
1270 : : !! mpi_enreg=information about MPI parallelization
1271 : : !! nfft=(effective) number of FFT grid points (for this processor)
1272 : : !! ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft
1273 : : !! frin(nfft)=real space input function
1274 : : !!
1275 : : !! OUTPUT
1276 : : !! frredgr(nfft,3)= reduced gradient of input function (same units as frin)
1277 : : !!
1278 : : !! SOURCE
1279 : :
1280 : 5 : subroutine redgr(frin,frredgr,mpi_enreg,nfft,ngfft)
1281 : :
1282 : : !Arguments ------------------------------------
1283 : : !scalars
1284 : : integer,intent(in) :: nfft
1285 : : type(MPI_type),intent(in) :: mpi_enreg
1286 : : !arrays
1287 : : integer,intent(in) :: ngfft(18)
1288 : : real(dp),intent(in) :: frin(nfft)
1289 : : real(dp),intent(out) :: frredgr(nfft,3)
1290 : :
1291 : : !Local variables-------------------------------
1292 : : !scalars
1293 : : integer :: cplex_tmp,i1,i2,i3,id,idir,ifft,ig,ii,ing,n1,n2,n3
1294 : : !arrays
1295 : 5 : real(dp),allocatable :: gg(:,:),wkcmpx(:,:),work(:),workgr(:,:)
1296 : : ! *************************************************************************
1297 : :
1298 : : !Only real arrays are treated
1299 : 5 : cplex_tmp=1
1300 : :
1301 : : !Keep local copy of fft dimensions
1302 : 5 : n1=ngfft(1) ; n2=ngfft(2) ; n3=ngfft(3)
1303 : :
1304 : : !In order to speed the routine, precompute the components of g, including 2pi factor
1305 [ - + + - : 15 : ABI_MALLOC(gg,(max(n1,n2,n3),3))
+ - + - -
+ + - ]
1306 [ + + ]: 20 : do ii=1,3
1307 : 15 : id=ngfft(ii)/2+2
1308 [ + + ]: 263 : do ing=1,ngfft(ii)
1309 : 248 : ig=ing-(ing/id)*ngfft(ii)-1
1310 : 263 : gg(ing,ii)=two_pi*ig
1311 : : end do
1312 : : ! Note that the G <-> -G symmetry must be maintained
1313 [ + + ]: 20 : if(mod(ngfft(ii),2)==0)gg(ngfft(ii)/2+1,ii)=zero
1314 : : end do
1315 : :
1316 [ - + + - : 15 : ABI_MALLOC(wkcmpx,(2,nfft))
+ - + - -
+ - + ]
1317 [ - + + - : 15 : ABI_MALLOC(work,(nfft))
+ - + - -
+ - + ]
1318 [ - + + - : 10 : ABI_MALLOC(workgr,(2,nfft))
- + + - ]
1319 : :
1320 : : !Obtain rho(G) in wkcmpx from input rho(r)
1321 [ + + ]: 24544 : work(:)=frin(:)
1322 : :
1323 : 5 : call fourdp(cplex_tmp,wkcmpx,work,-1,mpi_enreg,nfft,1,ngfft,0)
1324 : :
1325 : : !Gradient calculation for three reduced components in turn.
1326 : : !Code duplicated to remove logic from loops.
1327 [ + + ]: 20 : do idir=1,3
1328 [ + + ]: 15 : if(idir==1) then
1329 : : !$OMP PARALLEL DO PRIVATE(ifft)
1330 [ + + ]: 91 : do i3=1,n3
1331 : 86 : ifft=(i3-1)*n1*n2
1332 [ + + ]: 1514 : do i2=1,n2
1333 [ + + ]: 26048 : do i1=1,n1
1334 : 24539 : ifft=ifft+1
1335 : : ! Multiply by i 2pi G(idir)
1336 : 24539 : workgr(2,ifft)= gg(i1,idir)*wkcmpx(1,ifft)
1337 : 25962 : workgr(1,ifft)=-gg(i1,idir)*wkcmpx(2,ifft)
1338 : : end do
1339 : : end do
1340 : : end do
1341 [ + + ]: 10 : else if(idir==2) then
1342 : : !$OMP PARALLEL DO PRIVATE(ifft)
1343 [ + + ]: 91 : do i3=1,n3
1344 : 86 : ifft=(i3-1)*n1*n2
1345 [ + + ]: 1514 : do i2=1,n2
1346 [ + + ]: 26048 : do i1=1,n1
1347 : 24539 : ifft=ifft+1
1348 : : ! Multiply by i 2pi G(idir)
1349 : 24539 : workgr(2,ifft)= gg(i2,idir)*wkcmpx(1,ifft)
1350 : 25962 : workgr(1,ifft)=-gg(i2,idir)*wkcmpx(2,ifft)
1351 : : end do
1352 : : end do
1353 : : end do
1354 : : else
1355 : : !$OMP PARALLEL DO PRIVATE(ifft)
1356 [ + + ]: 91 : do i3=1,n3
1357 : 86 : ifft=(i3-1)*n1*n2
1358 [ + + ]: 1514 : do i2=1,n2
1359 [ + + ]: 26048 : do i1=1,n1
1360 : 24539 : ifft=ifft+1
1361 : : ! Multiply by i 2pi G(idir)
1362 : 24539 : workgr(2,ifft)= gg(i3,idir)*wkcmpx(1,ifft)
1363 : 25962 : workgr(1,ifft)=-gg(i3,idir)*wkcmpx(2,ifft)
1364 : : end do
1365 : : end do
1366 : : end do
1367 : : end if !idir
1368 : :
1369 : 15 : call fourdp(cplex_tmp,workgr,work,1,mpi_enreg,nfft,1,ngfft,0)
1370 : :
1371 : : !$OMP PARALLEL DO
1372 [ + + ]: 73637 : do ifft=1,nfft
1373 : 73632 : frredgr(ifft,idir)=work(ifft)
1374 : : end do
1375 : :
1376 : : end do !idir
1377 : :
1378 [ - + ]: 5 : ABI_FREE(gg)
1379 [ - + ]: 5 : ABI_FREE(wkcmpx)
1380 [ - + ]: 5 : ABI_FREE(work)
1381 [ - + ]: 5 : ABI_FREE(workgr)
1382 : :
1383 : 5 : end subroutine redgr
1384 : : !!***
1385 : :
1386 : : !!****f* m_spacepar/hartrestr
1387 : : !! NAME
1388 : : !! hartrestr
1389 : : !!
1390 : : !! FUNCTION
1391 : : !! To be called for strain perturbation only
1392 : : !! Compute the inhomogenous terms generated by the strain derivative of
1393 : : !! Hartree potential due to the ground state charge rho(G)
1394 : : !!
1395 : : !! FFT of (rho(G)/pi)*[d(1/G**2)/d(strain) - delta(diagonal strain)*(1/G**2)]
1396 : : !!
1397 : : !! NOTES
1398 : : !! *based largely on hartre.f
1399 : : !! *Modified code to avoid if statements inside loops to skip G=0.
1400 : : !! Replaced if statement on G^2>gsqcut to skip G s outside where
1401 : : !! rho(G) should be 0. Effect is negligible but gsqcut should be
1402 : : !! used to be strictly consistent with usage elsewhere in code.
1403 : : !! *The speed-up is provided by doing a few precomputations outside
1404 : : !! the inner loop. One variable size array is needed for this (gq).
1405 : : !!
1406 : : !! INPUTS
1407 : : !! gsqcut=cutoff value on G**2 for sphere inside fft box.
1408 : : !! idir=direction of the current perturbation
1409 : : !! ipert=type of the perturbation
1410 : : !! mpi_enreg=information about MPI parallelization
1411 : : !! natom=number of atoms in cell.
1412 : : !! nfft=number of fft grid points (gsqcut=(boxcut**2)*ecut/(2._dp*(Pi**2))
1413 : : !! ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft
1414 : : !! rhog(2,nfft)=array for Fourier transform of GS electron density
1415 : : !! rprimd(3,3)=dimensional primitive translations in real space (bohr)
1416 : : !!
1417 : : !! OUTPUT
1418 : : !! vhartr1(nfft)=Inhomogeneous term in strain-perturbation-induced Hartree
1419 : : !! potential in real space,
1420 : : !!
1421 : : !! SOURCE
1422 : :
1423 : 2875 : subroutine hartrestr(gsqcut,idir,ipert,mpi_enreg,natom,nfft,ngfft,rhog,rprimd,vhartr1)
1424 : :
1425 : : !Arguments ------------------------------------
1426 : : !scalars
1427 : : integer,intent(in) :: idir,ipert,natom,nfft
1428 : : real(dp),intent(in) :: gsqcut
1429 : : type(MPI_type),intent(in) :: mpi_enreg
1430 : : !arrays
1431 : : integer,intent(in) :: ngfft(18)
1432 : : real(dp),intent(in) :: rhog(2,nfft),rprimd(3,3)
1433 : : real(dp),intent(out) :: vhartr1(nfft)
1434 : :
1435 : : !Local variables-------------------------------
1436 : : !scalars
1437 : : integer,parameter :: im=2,re=1
1438 : : integer :: i1,i2,i23,i3,id2,id3,ig,ig2,ig3,ii,ii1,ing,istr,ka,kb,n1,n2,n3
1439 : : real(dp),parameter :: tolfix=1.000000001_dp
1440 : : real(dp) :: cutoff,ddends,den,dgsds,gqg2p3,gqgm12,gqgm13,gqgm23,gs,gs2,gs3
1441 : : real(dp) :: term,ucvol
1442 : : character(len=500) :: message
1443 : : !arrays
1444 : : integer,save :: idx(12)=(/1,1,2,2,3,3,3,2,3,1,2,1/)
1445 : : integer :: id(3)
1446 : 2875 : integer, contiguous, pointer :: fftn2_distrib(:),ffti2_local(:)
1447 : 2875 : integer, contiguous, pointer :: fftn3_distrib(:),ffti3_local(:)
1448 : : real(dp) :: dgmetds(3,3),gmet(3,3),gprimd(3,3),gqr(3),rmet(3,3)
1449 : 2875 : real(dp),allocatable :: gq(:,:),work1(:,:)
1450 : : ! *************************************************************************
1451 : :
1452 [ - + ]: 1387 : if( .not. (ipert==natom+3 .or. ipert==natom+4))then
1453 : : write(message, '(a,i0,a,a)' )&
1454 : 0 : 'From the calling routine, ipert=',ipert,ch10,&
1455 : 0 : 'so this routine for the strain perturbation should not be called.'
1456 : 0 : ABI_BUG(message)
1457 : : end if
1458 : :
1459 : 2875 : call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
1460 : :
1461 : 2875 : n1=ngfft(1) ; n2=ngfft(2) ; n3=ngfft(3)
1462 : :
1463 : : !Get the distrib associated with this fft_grid
1464 : 2875 : call ptabs_fourdp(mpi_enreg,n2,n3,fftn2_distrib,ffti2_local,fftn3_distrib,ffti3_local)
1465 : :
1466 : : !Initialize a few quantities
1467 : 2875 : cutoff=gsqcut*tolfix
1468 : :
1469 : 2875 : istr=idir + 3*(ipert-natom-3)
1470 : :
1471 [ - + ]: 2875 : if(istr<1 .or. istr>6)then
1472 : : write(message, '(a,i10,a,a,a)' )&
1473 : 0 : & 'Input dir gives istr=',istr,' not allowed.',ch10,&
1474 : 0 : & 'Possible values are 1,2,3,4,5,6 only.'
1475 : 0 : ABI_BUG(message)
1476 : : end if
1477 : :
1478 : 2875 : ka=idx(2*istr-1);kb=idx(2*istr)
1479 [ + + ]: 11500 : do ii = 1,3
1480 [ + + ]: 37375 : dgmetds(:,ii)=-(gprimd(ka,:)*gprimd(kb,ii)+gprimd(kb,:)*gprimd(ka,ii))
1481 : : end do
1482 : : !For historical reasons:
1483 [ + + + + ]: 37375 : dgmetds(:,:)=0.5_dp*dgmetds(:,:)
1484 : :
1485 : : !In order to speed the routine, precompute the components of g+q
1486 : : !Also check if the booked space was large enough...
1487 [ - + + - : 8625 : ABI_MALLOC(gq,(3,max(n1,n2,n3)))
+ - + - -
+ + - ]
1488 [ + + ]: 11500 : do ii=1,3
1489 : 8625 : id(ii)=ngfft(ii)/2+2
1490 [ + + ]: 170946 : do ing=1,ngfft(ii)
1491 : 159446 : ig=ing-(ing/id(ii))*ngfft(ii)-1
1492 : 168071 : gq(ii,ing)=ig
1493 : : end do
1494 : : end do
1495 : :
1496 [ - + + - : 8625 : ABI_MALLOC(work1,(2,nfft))
+ - + - -
+ - + ]
1497 : 2875 : id2=n2/2+2
1498 : 2875 : id3=n3/2+2
1499 : : !Triple loop on each dimension
1500 [ + + ]: 67797 : do i3=1,n3
1501 : 64922 : ig3=i3-(i3/id3)*n3-1
1502 : : ! Precompute some products that do not depend on i2 and i1
1503 : 64922 : gqr(3)=gq(3,i3)
1504 : 64922 : gs3=gq(3,i3)*gq(3,i3)*gmet(3,3)
1505 : 64922 : gqgm23=gq(3,i3)*gmet(2,3)*2
1506 : 64922 : gqgm13=gq(3,i3)*gmet(1,3)*2
1507 : :
1508 [ + + ]: 1177463 : do i2=1,n2
1509 [ + - ]: 1174588 : if (fftn2_distrib(i2)==mpi_enreg%me_fft) then
1510 : 1109666 : gqr(2)=gq(2,i2)
1511 : 1109666 : gs2=gs3+ gq(2,i2)*(gq(2,i2)*gmet(2,2)+gqgm23)
1512 : 1109666 : gqgm12=gq(2,i2)*gmet(1,2)*2
1513 : 1109666 : gqg2p3=gqgm13+gqgm12
1514 : 1109666 : ig2=i2-(i2/id2)*n2-1
1515 : : ! i23=n1*((i2-1)+n2*(i3-1))
1516 : 1109666 : i23=n1*((ffti2_local(i2)-1)+(n2/mpi_enreg%nproc_fft)*(i3-1))
1517 : : ! Do the test that eliminates the Gamma point outside
1518 : : ! of the inner loop
1519 : 1109666 : ii1=1
1520 [ + + ]: 1109666 : if(i23==0 .and. ig2==0 .and. ig3==0)then
1521 : 2875 : ii1=2
1522 : 2875 : work1(re,1+i23)=0.0_dp
1523 : 2875 : work1(im,1+i23)=0.0_dp
1524 : : end if
1525 : :
1526 : : ! Final inner loop on the first dimension
1527 : : ! (note the lower limit)
1528 [ + + ]: 20560083 : do i1=ii1,n1
1529 : 19450417 : gs=gs2+ gq(1,i1)*(gq(1,i1)*gmet(1,1)+gqg2p3)
1530 : 19450417 : ii=i1+i23
1531 [ + + ]: 20560083 : if(gs<=cutoff)then
1532 : 6455590 : den=piinv/gs
1533 : 6455590 : gqr(1)=gq(1,i1)
1534 : : dgsds=&
1535 : : & (gqr(1)*(dgmetds(1,1)*gqr(1)+dgmetds(1,2)*gqr(2)+dgmetds(1,3)*gqr(3))+ &
1536 : : & gqr(2)*(dgmetds(2,1)*gqr(1)+dgmetds(2,2)*gqr(2)+dgmetds(2,3)*gqr(3))+ &
1537 : 6455590 : & gqr(3)*(dgmetds(3,1)*gqr(1)+dgmetds(3,2)*gqr(2)+dgmetds(3,3)*gqr(3)) )
1538 : 6455590 : ddends=-piinv*dgsds/gs**2
1539 [ + + ]: 6455590 : if(istr<=3)then
1540 : 3397998 : term=2.0_dp*ddends-den
1541 : : else
1542 : 3057592 : term=2.0_dp*ddends
1543 : : end if
1544 : 6455590 : work1(re,ii)=rhog(re,ii)*term
1545 : 6455590 : work1(im,ii)=rhog(im,ii)*term
1546 : : else
1547 : 12994827 : work1(re,ii)=0.0_dp
1548 : 12994827 : work1(im,ii)=0.0_dp
1549 : : end if
1550 : :
1551 : : end do ! End loop on i1
1552 : : end if
1553 : : end do ! End loop on i2
1554 : : end do ! End loop on i3
1555 : :
1556 [ - + ]: 2875 : ABI_FREE(gq)
1557 : :
1558 : : !Fourier Transform Vhartree.
1559 : : !Vh in reciprocal space was stored in work1
1560 : 2875 : call fourdp(1,work1,vhartr1,1,mpi_enreg,nfft,1,ngfft,0)
1561 : :
1562 [ - + ]: 2875 : ABI_FREE(work1)
1563 : :
1564 [ + + ]: 2875 : end subroutine hartrestr
1565 : : !!***
1566 : :
1567 : : !!****f* m_spacepar/symrhg
1568 : : !! NAME
1569 : : !! symrhg
1570 : : !!
1571 : : !! FUNCTION
1572 : : !! From rho(r), generate rho(G), symmetrize it, and
1573 : : !! come back to the real space for a symmetrized rho(r).
1574 : : !!
1575 : : !! INPUTS
1576 : : !! cplex=1 if rhor is real, 2 if rhor is complex
1577 : : !! gprimd(3,3)=dimensional reciprocal space primitive translations
1578 : : !! irrzon(nfft,2,(nspden/nsppol)-3*(nspden/4))=irreducible zone data
1579 : : !! mpi_enreg=information about MPI parallelization
1580 : : !! nfft=(effective) number of FFT grid points (for this processor)
1581 : : !! ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft
1582 : : !! nspden=number of spin-density components
1583 : : !! nsppol=1 for unpolarized, 2 for spin-polarized
1584 : : !! nsym=number of symmetry elements.
1585 : : !! phnons(2,nfft,(nspden/nsppol)-3*(nspden/4))=nonsymmorphic translation phases
1586 : : !! rprimd(3,3)=dimensional real space primitive translations
1587 : : !! symafm(nsym)=(anti)ferromagnetic part of symmetry operations
1588 : : !! symrel(3,3,nsym)=symmetry matrices in real space (integers)
1589 : : !! tnons(3,nsym)=reduced nonsymmorphic translations
1590 : : !!
1591 : : !! OUTPUT
1592 : : !! rhog(2,nfft)=symmetrized rho(G) (total) electron density in G space
1593 : : !!
1594 : : !! SIDE EFFECTS
1595 : : !! Input/Output
1596 : : !! rhor(cplex*nfft,nspden)=array for electron density in electrons/bohr**3.
1597 : : !! Input, but also output, if symmetrization is applied.
1598 : : !! Also output if nspden > 1 (change spin components)
1599 : : !!
1600 : : !! NOTES
1601 : : !! When using spin-polarization (nspden==2),
1602 : : !! put total density in first half of rhor array and spin up in second half
1603 : : !! If (nspden=2 and nsppol=2) the density is transformed as (up,down) => (up+down,up)
1604 : : !! If (nspden=2 and nsppol=1) anti-ferromagnetic symmetry operations
1605 : : !! must be used, such as to transform (2*up) => (up+down,up)
1606 : : !! In spin-polarized, and if there is no symmetry to be
1607 : : !! applied on the system, only the total density is generated in G space
1608 : : !!
1609 : : !! SOURCE
1610 : :
1611 : 72047 : subroutine symrhg(cplex,gprimd,irrzon,mpi_enreg,nfft,nfftot,ngfft,nspden,nsppol,nsym,&
1612 : 72047 : & phnons,rhog,rhor,rprimd,symafm,symrel,tnons)
1613 : :
1614 : : !Arguments ------------------------------------
1615 : : !scalars
1616 : : integer,intent(in) :: cplex,nfft,nfftot,nspden,nsppol,nsym
1617 : : type(MPI_type),intent(in) :: mpi_enreg
1618 : : !arrays
1619 : : integer,intent(in) :: irrzon(nfftot**(1-1/nsym),2,(nspden/nsppol)-3*(nspden/4)),ngfft(18)
1620 : : integer,intent(in) :: symafm(nsym),symrel(3,3,nsym)
1621 : : real(dp),intent(in) :: gprimd(3,3),phnons(2,nfftot**(1-1/nsym),(nspden/nsppol)-3*(nspden/4)),rprimd(3,3)
1622 : : real(dp),intent(inout) :: rhor(cplex*nfft,nspden)
1623 : : real(dp),intent(out) :: rhog(2,nfft)
1624 : : real(dp),intent(in) :: tnons(3,nsym)
1625 : :
1626 : : !Local variables-------------------------------
1627 : : !scalars
1628 : : integer :: id1,id2,id3,ier,imagn,ind,ind2,indsy,ispden,isym,iup,izone,izone_max,j,j1,j2,j3,jsym
1629 : : integer :: k1,k2,k3,l1,l2,l3,me_fft
1630 : : integer :: n1,n2,n3,nd2,nproc_fft,nspden_eff,nsym_used,numpt,nup
1631 : : integer :: r2,rep,spaceComm
1632 : : logical,parameter :: afm_noncoll=.true. ! TRUE if antiferro symmetries are used in non-collinear magnetism
1633 : : real(dp) :: arg,tau1,tau2,tau3
1634 : : real(dp) :: magxsu1,magxsu2,magysu1,magysu2,magzsu1,magzsu2,mxi,mxr,myi,myr,mzi,mzr,phi,phr,rhosu1,rhosu2
1635 : : !character(len=500) :: message
1636 : : !arrays
1637 : 72047 : integer,allocatable :: isymg(:)
1638 : 72047 : integer, contiguous, pointer :: fftn2_distrib(:),ffti2_local(:)
1639 : 72047 : integer, contiguous, pointer :: fftn3_distrib(:),ffti3_local(:)
1640 : : real(dp) :: tsec(2)
1641 : 72047 : real(dp),allocatable :: magngx(:,:),magngy(:,:),magngz(:,:)
1642 : 72047 : real(dp),allocatable :: rhosu1_arr(:),rhosu2_arr(:),work(:)
1643 : 72047 : real(dp),allocatable :: symafm_used(:),symrec_cart(:,:,:),symrel_cart(:,:,:),tnons_used(:,:),sym_det(:)
1644 : : !*************************************************************************
1645 : : !
1646 : : !Note the timing channel 17 excludes the different Fourier transforms
1647 : :
1648 [ + - + - : 216141 : ABI_MALLOC(work,(cplex*nfft))
+ - - + -
+ ]
1649 : :
1650 : : !Special treatment for spin-polarized case
1651 [ + + + + ]: 72047 : if(nspden==2 .and. nsppol==2) then
1652 : : ! When nspden=2 and nsppol=2, put total density in first half
1653 : : ! of rhor array and spin up in second half (up,down) => (up+down,up)
1654 : 6610 : call timab(17,1,tsec)
1655 [ + + ]: 115874757 : work(:)=rhor(:,1) ! up => work
1656 [ + + ]: 115874757 : rhor(:,1)=rhor(:,1)+rhor(:,2) ! up+down
1657 [ + + ]: 115874757 : rhor(:,2)=work(:) ! work => up
1658 : 6610 : call timab(17,2,tsec)
1659 : : end if
1660 : :
1661 : : !Special treatment for antiferromagnetism case
1662 [ + + + + ]: 72047 : if(nspden==2 .and. nsppol==1) then
1663 : 423 : call timab(17,1,tsec)
1664 : : ! When nspden=2 and nsppol=1, (2*up) => (2*up,up)
1665 : : ! Indeed, what was delivered to the present routine is a "total" density,
1666 : : ! obtained from occupation numbers varying between 0 and 2,
1667 : : ! but for spin up only potential.
1668 [ + + ]: 6079543 : rhor(:,2)=half*rhor(:,1)
1669 : 423 : call timab(17,2,tsec)
1670 : : end if
1671 : :
1672 : : !Special treatment for non-collinear magnetism case
1673 [ + + ]: 72047 : if(nspden==4) then
1674 : 2872 : call timab(17,1,tsec)
1675 : : !FR the half factors missing are recovered in dfpt_mkvxc_noncoll and dfpt_accrho
1676 [ + + ]: 26180531 : rhor(:,1)=rhor(:,1)+rhor(:,4) !nup+ndown
1677 [ + + ]: 26180531 : rhor(:,2)=rhor(:,2)-rhor(:,1) !mx (n+mx-n)
1678 [ + + ]: 26180531 : rhor(:,3)=rhor(:,3)-rhor(:,1) !my (n+my-n)
1679 [ + + ]: 26180531 : rhor(:,4)=rhor(:,1)-two*rhor(:,4) !mz=n-2ndown
1680 : 2872 : call timab(17,2,tsec)
1681 : : end if
1682 : :
1683 : :
1684 [ + + ]: 72047 : if(nsym==1)then
1685 : :
1686 [ + + - + ]: 26230 : if(nspden==2 .and. nsppol==1) then ! There must be at least one anti-ferromagnetic operation
1687 : 0 : ABI_BUG('In the antiferromagnetic case, nsym cannot be 1')
1688 : : end if
1689 : :
1690 : : ! If not using symmetry, still want total density in G space rho(G).
1691 : : ! Fourier transform (incl normalization) to get rho(G)
1692 [ + + ]: 369568944 : work(:)=rhor(:,1)
1693 : 26230 : call fourdp(cplex,rhog,work,-1,mpi_enreg,nfft,1,ngfft,0)
1694 : : else
1695 : :
1696 : : ! Treat either full density, spin-up density or magnetization
1697 : : ! Note the decrease of ispden to the value 1, in order to finish
1698 : : ! with rhog of the total density (and not the spin-up density or magnetization)
1699 [ + + ]: 45817 : nspden_eff=nspden;if (nspden==4) nspden_eff=1
1700 [ + + ]: 96484 : do ispden=nspden_eff,1,-1
1701 : :
1702 : : ! Prepare the density to be symmetrized, in the reciprocal space
1703 [ + + + + ]: 50667 : if(nspden==1 .or. nsppol==2 .or. (nspden==4.and.(.not.afm_noncoll)))then
1704 : 1108436 : imagn=1
1705 : : nsym_used=0
1706 [ + + ]: 1108436 : do isym=1,nsym
1707 [ + + ]: 1108436 : if(symafm(isym)==1)nsym_used=nsym_used+1
1708 : : ! DEBUG
1709 : : ! write(std_out,*)' symrhg : isym,symafm(isym)',isym,symafm(isym)
1710 : : ! ENDDEBUG
1711 : : end do
1712 [ + + + - ]: 2657 : else if(nspden==2 .and. nsppol==1)then ! antiferromagnetic case
1713 : 846 : imagn=ispden
1714 : 846 : nsym_used=nsym/ispden
1715 [ + - ]: 1811 : else if (nspden==4) then
1716 : 1811 : imagn=1
1717 : 1811 : nsym_used=nsym/ispden
1718 : : end if
1719 : :
1720 : : ! write(std_out,*)' symrhg : nsym_used=',nsym_used
1721 : :
1722 : : ! rhor -fft-> rhog (rhog is used as work space)
1723 : : ! Note : it should be possible to reuse rhog in the antiferromagnetic case this would avoid one FFT
1724 [ + + ]: 755152265 : work(:)=rhor(:,ispden)
1725 : 50667 : call fourdp(cplex,rhog,work,-1,mpi_enreg,nfft,1,ngfft,0)
1726 [ + + ]: 50667 : if (nspden==4) then
1727 [ - + + - : 5433 : ABI_MALLOC(magngx,(2,nfft))
+ - + - -
+ - + -
+ ]
1728 [ - + + - : 3622 : ABI_MALLOC(magngy,(2,nfft))
- + - + -
+ ]
1729 [ - + + - : 3622 : ABI_MALLOC(magngz,(2,nfft))
- + - + -
+ ]
1730 [ + + ]: 13871209 : work(:)=rhor(:,2)
1731 : 1811 : call fourdp(cplex,magngx,work,-1,mpi_enreg,nfft,1,ngfft,0)
1732 [ + + ]: 13871209 : work(:)=rhor(:,3)
1733 : 1811 : call fourdp(cplex,magngy,work,-1,mpi_enreg,nfft,1,ngfft,0)
1734 [ + + ]: 13871209 : work(:)=rhor(:,4)
1735 : 1811 : call fourdp(cplex,magngz,work,-1,mpi_enreg,nfft,1,ngfft,0)
1736 : : end if
1737 : :
1738 : : ! Begins the timing here only , to exclude FFTs
1739 : 50667 : call timab(17,1,tsec)
1740 : :
1741 : 50667 : n1=ngfft(1);n2=ngfft(2);n3=ngfft(3);nproc_fft=ngfft(10);me_fft=ngfft(11);nd2=n2/nproc_fft
1742 : :
1743 : : ! Get the distrib associated with this fft_grid
1744 : 50667 : call ptabs_fourdp(mpi_enreg,n2,n3,fftn2_distrib,ffti2_local,fftn3_distrib,ffti3_local)
1745 : :
1746 : : ! The following is only valid for total, up or dn density
1747 : : ! -------------------------------------------------------
1748 : :
1749 : : ! Get maxvalue of izone
1750 [ + + + + ]: 760617440 : izone_max=count(irrzon(:,2,imagn)>0)
1751 [ + + + - : 152001 : ABI_MALLOC(rhosu1_arr,(izone_max))
+ - + + -
+ - + -
+ ]
1752 [ + + + - : 101334 : ABI_MALLOC(rhosu2_arr,(izone_max))
- + - + -
+ ]
1753 : :
1754 : 50667 : numpt=0
1755 [ + - ]: 111317996 : do izone=1,nfftot
1756 : :
1757 : : ! Get repetition number
1758 : 111317996 : rep=irrzon(izone,2,imagn)
1759 [ + + ]: 111317996 : if(rep==0)exit
1760 : :
1761 : : ! Compute number of unique points in this symm class:
1762 : 111267329 : nup=nsym_used/rep
1763 : :
1764 : : ! Accumulate charge over equivalent points
1765 : 111267329 : rhosu1=zero
1766 : 111267329 : rhosu2=zero
1767 [ + + ]: 871653951 : do iup=1,nup
1768 : 760386622 : ind=irrzon(iup+numpt,1,imagn)
1769 : 760386622 : j=ind-1;j1=modulo(j,n1);j2=modulo(j/n1,n2);j3=j/(n1*n2);
1770 [ + + ]: 871653951 : if(fftn2_distrib(j2+1)==me_fft) then ! this ind is to be treated by me_fft
1771 : 744440814 : r2=ffti2_local(j2+1) - 1
1772 : 744440814 : ind=n1*(nd2*j3+r2)+j1+1 !this is ind in the current proc
1773 : : rhosu1=rhosu1+rhog(1,ind)*phnons(1,iup+numpt,imagn)&
1774 : 744440814 : & -rhog(2,ind)*phnons(2,iup+numpt,imagn)
1775 : : rhosu2=rhosu2+rhog(2,ind)*phnons(1,iup+numpt,imagn)&
1776 : 744440814 : & +rhog(1,ind)*phnons(2,iup+numpt,imagn)
1777 : : end if
1778 : :
1779 : : end do
1780 : 111267329 : rhosu1=rhosu1/dble(nup)
1781 : 111267329 : rhosu2=rhosu2/dble(nup)
1782 : 111267329 : rhosu1_arr(izone)=rhosu1
1783 : 111267329 : rhosu2_arr(izone)=rhosu2
1784 : : ! Keep index of how many points have been considered:
1785 : 111317996 : numpt=numpt+nup
1786 : :
1787 : : end do ! End loop over izone
1788 : :
1789 : : ! Reduction in case of FFT parallelization
1790 [ + + ]: 50667 : if(mpi_enreg%nproc_fft>1)then
1791 : 2648 : spaceComm=mpi_enreg%comm_fft
1792 : 2648 : call xmpi_sum(rhosu1_arr,spaceComm,ier)
1793 : 2648 : call xmpi_sum(rhosu2_arr,spaceComm,ier)
1794 : : end if
1795 : :
1796 : : ! Now symmetrize the density
1797 : : numpt=0
1798 [ + - ]: 111317996 : do izone=1,nfftot
1799 : :
1800 : : ! Get repetition number
1801 : 111317996 : rep=irrzon(izone,2,imagn)
1802 [ + + ]: 111317996 : if(rep==0)exit
1803 : :
1804 : : ! Compute number of unique points in this symm class:
1805 : 111267329 : nup=nsym_used/rep
1806 : :
1807 : : ! Define symmetrized rho(G) at equivalent points:
1808 [ + + ]: 871653951 : do iup=1,nup
1809 : 760386622 : ind=irrzon(iup+numpt,1,imagn)
1810 : : ! decompose ind-1=n1(n2 j3+ j2)+j1
1811 : 760386622 : j=ind-1;j1=modulo(j,n1);j2=modulo(j/n1,n2);j3=j/(n1*n2);
1812 [ + + ]: 871653951 : if(fftn2_distrib(j2+1)==me_fft) then ! this ind is to be treated by me_fft
1813 : 744440814 : r2=ffti2_local(j2+1) - 1
1814 : : ! ind in the proc ind-1=n1(nd2 j3+ r2)+j1
1815 : 744440814 : ind=n1*(nd2*j3+r2)+j1+1 !this is ind in the current proc
1816 : : rhog(1,ind)=rhosu1_arr(izone)*phnons(1,iup+numpt,imagn)&
1817 : 744440814 : & +rhosu2_arr(izone)*phnons(2,iup+numpt,imagn)
1818 : : rhog(2,ind)=rhosu2_arr(izone)*phnons(1,iup+numpt,imagn)&
1819 : 744440814 : & -rhosu1_arr(izone)*phnons(2,iup+numpt,imagn)
1820 : : end if
1821 : : end do
1822 : :
1823 : : ! Keep index of how many points have been considered:
1824 : 111317996 : numpt=numpt+nup
1825 : :
1826 : : end do ! End loop over izone
1827 : :
1828 [ - + ]: 50667 : ABI_FREE(rhosu1_arr)
1829 [ - + ]: 50667 : ABI_FREE(rhosu2_arr)
1830 : :
1831 : : ! The following is only valid for magnetization
1832 : : ! ---------------------------------------------
1833 [ + + ]: 50667 : if (nspden==4) then
1834 : :
1835 : 1811 : id1=n1/2+2
1836 : 1811 : id2=n2/2+2
1837 : 1811 : id3=n3/2+2
1838 : :
1839 : : ! Transfer symmetries in cartesian coordinates
1840 : : ! Compute symmetries in reciprocal space in cartesian coordinates
1841 [ - + + - : 5433 : ABI_MALLOC(symrec_cart,(3,3,nsym_used))
+ - + - -
+ - + ]
1842 [ - + + - : 3622 : ABI_MALLOC(symrel_cart,(3,3,nsym_used))
- + - + ]
1843 [ - + + - : 5433 : ABI_MALLOC(symafm_used,(nsym_used))
+ - + - -
+ - + ]
1844 [ - + + - : 5433 : ABI_MALLOC(tnons_used,(3,nsym_used))
+ - + - -
+ - + ]
1845 [ - + + - : 3622 : ABI_MALLOC(sym_det,(nsym_used))
- + - + ]
1846 : 20319 : jsym=0
1847 [ + + ]: 20319 : do isym=1,nsym
1848 : : if (symafm(isym)/=1.and.(.not.afm_noncoll)) cycle
1849 : 18508 : jsym=jsym+1
1850 [ + + ]: 74032 : tnons_used(:,jsym)=tnons(:,isym)
1851 : 18508 : symafm_used(jsym)=dble(symafm(isym))
1852 : 18508 : call symredcart(rprimd,gprimd,symrel_cart(:,:,jsym),symrel(:,:,isym))
1853 : 18508 : call matr3inv(symrel_cart(:,:,jsym),symrec_cart(:,:,jsym))
1854 : : sym_det(jsym) = symrel_cart(1,1,isym)*symrel_cart(2,2,isym)*symrel_cart(3,3,isym)+&
1855 : : & symrel_cart(2,1,isym)*symrel_cart(3,2,isym)*symrel_cart(1,3,isym)+&
1856 : : & symrel_cart(1,2,isym)*symrel_cart(2,3,isym)*symrel_cart(3,1,isym) - &
1857 : : & (symrel_cart(3,1,isym)*symrel_cart(2,2,isym)*symrel_cart(1,3,isym)+&
1858 : : & symrel_cart(2,1,isym)*symrel_cart(1,2,isym)*symrel_cart(3,3,isym)+&
1859 : 20319 : & symrel_cart(3,2,isym)*symrel_cart(2,3,isym)*symrel_cart(1,1,isym))
1860 : : end do
1861 : :
1862 [ + + + - ]: 13994089 : numpt=count(irrzon(:,1,imagn)>0)
1863 [ - + + - : 5433 : ABI_MALLOC(isymg,(numpt))
+ - + - -
+ - + ]
1864 [ + + ]: 13994089 : isymg=0
1865 [ + - - + : 5433 : ABI_MALLOC(rhosu1_arr,(3*izone_max))
+ - + - +
- - + -
+ ]
1866 [ - + + - : 3622 : ABI_MALLOC(rhosu2_arr,(3*izone_max))
- + - + ]
1867 : :
1868 : : ! Accumulate magnetization over equivalent points
1869 : : ! Use all symmetries (not only those linking different g points)
1870 : : ! Use Inverse[Transpose[symrel]]=symrec
1871 : 1811 : numpt=0
1872 [ + + ]: 4688121 : do izone=1,izone_max
1873 : 4686310 : magxsu1=zero;magxsu2=zero
1874 : 4686310 : magysu1=zero;magysu2=zero
1875 : 4686310 : magzsu1=zero;magzsu2=zero
1876 : 4686310 : ind=irrzon(1+numpt,1,1)
1877 : 4686310 : rep=irrzon(izone,2,1)
1878 : 4686310 : nup=nsym_used/rep
1879 : : ! Get coordinates in the range [0,n-1]
1880 : 4686310 : j=ind-1;l1=modulo(j,n1);l2=modulo(j/n1,n2);l3=j/(n1*n2)
1881 : : ! Get location of G vector (grid point) centered at 0 0 0
1882 : : !TO BE UNCOMMENTED
1883 : 4686310 : l3=l3-(l3/id3)*n3
1884 : 4686310 : l2=l2-(l2/id2)*n2
1885 : 4686310 : l1=l1-(l1/id1)*n1
1886 : :
1887 : 4686310 : jsym=0
1888 [ + + ]: 21116214 : do isym=1,nsym
1889 : : if (symafm(isym)/=1.and.(.not.afm_noncoll)) cycle
1890 : 16429904 : jsym=jsym+1
1891 : : ! The G vectors should transform as vectors in reciprocal space
1892 : : ! However, one acts with the INVERSE of the symmetry operation => Inverse[symrec]=Transpose[symrel]
1893 : 16429904 : j1=symrel(1,1,isym)*l1+symrel(2,1,isym)*l2+symrel(3,1,isym)*l3
1894 : 16429904 : j2=symrel(1,2,isym)*l1+symrel(2,2,isym)*l2+symrel(3,2,isym)*l3
1895 : 16429904 : j3=symrel(1,3,isym)*l1+symrel(2,3,isym)*l2+symrel(3,3,isym)*l3
1896 : 16429904 : k1=map_symrhg(j1,n1);k2=map_symrhg(j2,n2);k3=map_symrhg(j3,n3)
1897 : 16429904 : indsy=1+k1+n1*(k2+n2*k3)
1898 : 16429904 : ind2=-1;iup=numpt
1899 [ + + - + ]: 76629240 : do while (ind2/=indsy.and.iup<numpt+nup)
1900 : 60199336 : iup=iup+1;ind2=irrzon(iup,1,1)
1901 : : end do
1902 [ - + ]: 16429904 : if (ind2/=indsy) then
1903 : 0 : ABI_ERROR("ind2/=indsy in symrhg !")
1904 : : end if
1905 [ + + ]: 16429904 : if (isymg(iup)==0) isymg(iup)=jsym
1906 [ + + ]: 21116214 : if(fftn2_distrib(modulo((indsy-1)/n1,n2) + 1) == me_fft ) then ! this is indsy is to be treated by me_fft
1907 : 16253504 : indsy=n1*(nd2*k3+ ffti2_local(k2+1) -1)+k1+1 ! this is indsy in the current proc
1908 : :
1909 : : ! Working on this: the present coding will be detrimental for speed ! cos and sin are recomputed many times !
1910 : 16253504 : tau1=tnons_used(1,jsym)
1911 : 16253504 : tau2=tnons_used(2,jsym)
1912 : 16253504 : tau3=tnons_used(3,jsym)
1913 [ + + + - : 16253504 : if (abs(tau1)>tol12.or.abs(tau2)>tol12.or.abs(tau3)>tol12) then
+ + ]
1914 : : ! Compute exp(-2*Pi*I*G dot tau) using original G (equivalent of phnons in the collinear case)
1915 : 1288860 : arg=two_pi*(dble(l1)*tau1+dble(l2)*tau2+dble(l3)*tau3)
1916 : 1288860 : phr=cos(arg)
1917 : 1288860 : phi=-sin(arg)
1918 : : else
1919 : : phr=one
1920 : : phi=zero
1921 : : end if
1922 : 16253504 : phr=phr*symafm_used(jsym)
1923 : 16253504 : phi=phi*symafm_used(jsym)
1924 : : !TO BE COMMENTED
1925 : : ! phr=phnons(1,iup,imagn);if (rep==1) phr=phr*symafm_used(jsym) !if rep==2, symafm is already included in phnons
1926 : : ! phi=phnons(2,iup,imagn);if (rep==1) phi=phi*symafm_used(jsym) !(see irrzg.F90)
1927 : :
1928 : : ! The magnetization should transform as a vector in real space
1929 : : ! However, one acts with the INVERSE of the symmetry operation.
1930 : : ! => Inverse[symrel_cart] = Transpose[symrel_cart] because symrel_cart is unitary ?!?!?
1931 : 16253504 : mxr=sym_det(jsym)*(symrel_cart(1,1,jsym)*magngx(1,indsy)+symrel_cart(1,2,jsym)*magngy(1,indsy)+symrel_cart(1,3,jsym)*magngz(1,indsy))
1932 : 16253504 : mxi=sym_det(jsym)*(symrel_cart(1,1,jsym)*magngx(2,indsy)+symrel_cart(1,2,jsym)*magngy(2,indsy)+symrel_cart(1,3,jsym)*magngz(2,indsy))
1933 : 16253504 : myr=sym_det(jsym)*(symrel_cart(2,1,jsym)*magngx(1,indsy)+symrel_cart(2,2,jsym)*magngy(1,indsy)+symrel_cart(2,3,jsym)*magngz(1,indsy))
1934 : 16253504 : myi=sym_det(jsym)*(symrel_cart(2,1,jsym)*magngx(2,indsy)+symrel_cart(2,2,jsym)*magngy(2,indsy)+symrel_cart(2,3,jsym)*magngz(2,indsy))
1935 : 16253504 : mzr=sym_det(jsym)*(symrel_cart(3,1,jsym)*magngx(1,indsy)+symrel_cart(3,2,jsym)*magngy(1,indsy)+symrel_cart(3,3,jsym)*magngz(1,indsy))
1936 : 16253504 : mzi=sym_det(jsym)*(symrel_cart(3,1,jsym)*magngx(2,indsy)+symrel_cart(3,2,jsym)*magngy(2,indsy)+symrel_cart(3,3,jsym)*magngz(2,indsy))
1937 : :
1938 : : ! mxr=symrel_cart(1,1,jsym)*magngx(1,indsy)+symrel_cart(2,1,jsym)*magngy(1,indsy)+symrel_cart(3,1,jsym)*magngz(1,indsy)
1939 : : ! mxi=symrel_cart(1,1,jsym)*magngx(2,indsy)+symrel_cart(2,1,jsym)*magngy(2,indsy)+symrel_cart(3,1,jsym)*magngz(2,indsy)
1940 : : ! myr=symrel_cart(1,2,jsym)*magngx(1,indsy)+symrel_cart(2,2,jsym)*magngy(1,indsy)+symrel_cart(3,2,jsym)*magngz(1,indsy)
1941 : : ! myi=symrel_cart(1,2,jsym)*magngx(2,indsy)+symrel_cart(2,2,jsym)*magngy(2,indsy)+symrel_cart(3,2,jsym)*magngz(2,indsy)
1942 : : ! mzr=symrel_cart(1,3,jsym)*magngx(1,indsy)+symrel_cart(2,3,jsym)*magngy(1,indsy)+symrel_cart(3,3,jsym)*magngz(1,indsy)
1943 : : ! mzi=symrel_cart(1,3,jsym)*magngx(2,indsy)+symrel_cart(2,3,jsym)*magngy(2,indsy)+symrel_cart(3,3,jsym)*magngz(2,indsy)
1944 : :
1945 : 16253504 : magxsu1=magxsu1+mxr*phr-mxi*phi;magxsu2=magxsu2+mxi*phr+mxr*phi
1946 : 16253504 : magysu1=magysu1+myr*phr-myi*phi;magysu2=magysu2+myi*phr+myr*phi
1947 : 16253504 : magzsu1=magzsu1+mzr*phr-mzi*phi;magzsu2=magzsu2+mzi*phr+mzr*phi
1948 : : end if
1949 : : end do
1950 : 4686310 : rhosu1_arr(3*izone-2)=magxsu1/dble(nsym_used)
1951 : 4686310 : rhosu1_arr(3*izone-1)=magysu1/dble(nsym_used)
1952 : 4686310 : rhosu1_arr(3*izone )=magzsu1/dble(nsym_used)
1953 : 4686310 : rhosu2_arr(3*izone-2)=magxsu2/dble(nsym_used)
1954 : 4686310 : rhosu2_arr(3*izone-1)=magysu2/dble(nsym_used)
1955 : 4686310 : rhosu2_arr(3*izone )=magzsu2/dble(nsym_used)
1956 : 4688121 : numpt=numpt+nup
1957 : : end do
1958 : :
1959 : : ! Reduction in case of FFT parallelization
1960 [ + + ]: 1811 : if(mpi_enreg%nproc_fft>1)then
1961 : 60 : spaceComm=mpi_enreg%comm_fft
1962 : 60 : call xmpi_sum(rhosu1_arr,spaceComm,ier)
1963 : 60 : call xmpi_sum(rhosu2_arr,spaceComm,ier)
1964 : : end if
1965 : :
1966 : : ! Now symmetrize the magnetization at equivalent points
1967 : : ! Use Transpose[symrel]
1968 : : numpt=0
1969 [ + + ]: 4688121 : do izone=1,izone_max
1970 : 4686310 : rep=irrzon(izone,2,imagn)
1971 : 4686310 : nup=nsym_used/rep
1972 [ + + ]: 18678588 : do iup=1,nup
1973 : 13992278 : ind=irrzon(iup+numpt,1,imagn)
1974 : : ! Get coordinates in the range [0,n-1]
1975 : 13992278 : j=ind-1;j1=modulo(j,n1);j2=modulo(j/n1,n2);j3=j/(n1*n2)
1976 : : !TO BE UNCOMMENTED
1977 : : ! Get location of G vector (grid point) centered at 0 0 0
1978 : 13992278 : l3=j3-(j3/id3)*n3
1979 : 13992278 : l2=j2-(j2/id2)*n2
1980 : 13992278 : l1=j1-(j1/id1)*n1
1981 [ + + ]: 18678588 : if(fftn2_distrib(j2+1)==me_fft) then ! this ind is to be treated by me_fft
1982 : 13869398 : r2=ffti2_local(j2+1) - 1
1983 : 13869398 : ind=n1*(nd2*j3+r2)+j1+1 ! this is ind in the current proc
1984 : 13869398 : jsym=isymg(iup+numpt)
1985 [ - + ]: 13869398 : if (jsym==0) then
1986 : 0 : ABI_ERROR("jsym=0 in symrhg !")
1987 : : end if
1988 : 13869398 : magxsu1=rhosu1_arr(3*izone-2);magxsu2=rhosu2_arr(3*izone-2)
1989 : 13869398 : magysu1=rhosu1_arr(3*izone-1);magysu2=rhosu2_arr(3*izone-1)
1990 : 13869398 : magzsu1=rhosu1_arr(3*izone );magzsu2=rhosu2_arr(3*izone )
1991 : : ! Working on this: the present coding will be detrimental for speed ! cos and sin are recomputed many times !
1992 : 13869398 : tau1=tnons_used(1,jsym)
1993 : 13869398 : tau2=tnons_used(2,jsym)
1994 : 13869398 : tau3=tnons_used(3,jsym)
1995 [ + + + - : 13869398 : if (abs(tau1)>tol12.or.abs(tau2)>tol12.or.abs(tau3)>tol12) then
- + ]
1996 : : ! Compute exp(-2*Pi*I*G dot tau) using original G (equivalent of phnons in the collinear case)
1997 : 850426 : arg=two_pi*(dble(l1)*tau1+dble(l2)*tau2+dble(l3)*tau3)
1998 : 850426 : phr=cos(arg)
1999 : 850426 : phi=-sin(arg)
2000 : : else
2001 : : phr=one
2002 : : phi=zero
2003 : : end if
2004 : 13869398 : phr=phr*symafm_used(jsym)
2005 : 13869398 : phi=phi*symafm_used(jsym)
2006 : : !TO BE COMMENTED
2007 : : ! phr=phnons(1,iup,imagn);if (rep==1) phr=phr*symafm_used(jsym) !if rep==2, symafm is already included in phnons
2008 : : ! phi=phnons(2,iup,imagn);if (rep==1) phi=phi*symafm_used(jsym) !(see irrzg.F90)
2009 : : ! The magnetization should transform as a vector in real space
2010 : : ! => symrel_cart ?!?
2011 : 13869398 : mxr=sym_det(jsym)*(symrec_cart(1,1,jsym)*magxsu1+symrec_cart(2,1,jsym)*magysu1+symrec_cart(3,1,jsym)*magzsu1)
2012 : 13869398 : mxi=sym_det(jsym)*(symrec_cart(1,1,jsym)*magxsu2+symrec_cart(2,1,jsym)*magysu2+symrec_cart(3,1,jsym)*magzsu2)
2013 : 13869398 : myr=sym_det(jsym)*(symrec_cart(1,2,jsym)*magxsu1+symrec_cart(2,2,jsym)*magysu1+symrec_cart(3,2,jsym)*magzsu1)
2014 : 13869398 : myi=sym_det(jsym)*(symrec_cart(1,2,jsym)*magxsu2+symrec_cart(2,2,jsym)*magysu2+symrec_cart(3,2,jsym)*magzsu2)
2015 : 13869398 : mzr=sym_det(jsym)*(symrec_cart(1,3,jsym)*magxsu1+symrec_cart(2,3,jsym)*magysu1+symrec_cart(3,3,jsym)*magzsu1)
2016 : 13869398 : mzi=sym_det(jsym)*(symrec_cart(1,3,jsym)*magxsu2+symrec_cart(2,3,jsym)*magysu2+symrec_cart(3,3,jsym)*magzsu2)
2017 : : ! mxr=symrel_cart(1,1,jsym)*magxsu1+symrel_cart(1,2,jsym)*magysu1+symrel_cart(1,3,jsym)*magzsu1
2018 : : ! mxi=symrel_cart(1,1,jsym)*magxsu2+symrel_cart(1,2,jsym)*magysu2+symrel_cart(1,3,jsym)*magzsu2
2019 : : ! myr=symrel_cart(2,1,jsym)*magxsu1+symrel_cart(2,2,jsym)*magysu1+symrel_cart(2,3,jsym)*magzsu1
2020 : : ! myi=symrel_cart(2,1,jsym)*magxsu2+symrel_cart(2,2,jsym)*magysu2+symrel_cart(2,3,jsym)*magzsu2
2021 : : ! mzr=symrel_cart(3,1,jsym)*magxsu1+symrel_cart(3,2,jsym)*magysu1+symrel_cart(3,3,jsym)*magzsu1
2022 : : ! mzi=symrel_cart(3,1,jsym)*magxsu2+symrel_cart(3,2,jsym)*magysu2+symrel_cart(3,3,jsym)*magzsu2
2023 : 13869398 : magngx(1,ind)=mxr*phr-mxi*phi
2024 : 13869398 : magngx(2,ind)=mxi*phr+mxr*phi
2025 : 13869398 : magngy(1,ind)=myr*phr-myi*phi
2026 : 13869398 : magngy(2,ind)=myi*phr+myr*phi
2027 : 13869398 : magngz(1,ind)=mzr*phr-mzi*phi
2028 : 13869398 : magngz(2,ind)=mzi*phr+mzr*phi
2029 : : end if
2030 : : end do
2031 : 4688121 : numpt=numpt+nup
2032 : : end do
2033 [ - + ]: 1811 : ABI_FREE(isymg)
2034 [ - + ]: 1811 : ABI_FREE(rhosu1_arr)
2035 [ - + ]: 1811 : ABI_FREE(rhosu2_arr)
2036 [ - + ]: 1811 : ABI_FREE(symrec_cart)
2037 [ - + ]: 1811 : ABI_FREE(symrel_cart)
2038 [ - + ]: 1811 : ABI_FREE(sym_det)
2039 [ - + ]: 1811 : ABI_FREE(symafm_used)
2040 [ - + ]: 1811 : ABI_FREE(tnons_used)
2041 : :
2042 : : end if ! nspden==4
2043 : :
2044 : 50667 : call timab(17,2,tsec)
2045 : :
2046 : : ! Pull out full or spin up density, now symmetrized
2047 : 50667 : call fourdp(cplex,rhog,work,1,mpi_enreg,nfft,1,ngfft,0)
2048 [ + + ]: 755152265 : rhor(:,ispden)=work(:)
2049 [ + + ]: 96484 : if (nspden==4) then
2050 : 1811 : call fourdp(cplex,magngx,work,1,mpi_enreg,nfft,1,ngfft,0)
2051 [ + + ]: 13871209 : rhor(:,2)=work(:)
2052 : 1811 : call fourdp(cplex,magngy,work,1,mpi_enreg,nfft,1,ngfft,0)
2053 [ + + ]: 13871209 : rhor(:,3)=work(:)
2054 : 1811 : call fourdp(cplex,magngz,work,1,mpi_enreg,nfft,1,ngfft,0)
2055 [ + + ]: 13871209 : rhor(:,4)=work(:)
2056 [ - + ]: 1811 : ABI_FREE(magngx)
2057 [ - + ]: 1811 : ABI_FREE(magngy)
2058 [ - + ]: 1811 : ABI_FREE(magngz)
2059 : : end if
2060 : :
2061 : : end do ! ispden
2062 : :
2063 : : end if ! End on the condition nsym==1
2064 : :
2065 [ - + - + : 144094 : ABI_FREE(work)
- + - + -
+ - + -
+ ]
2066 : :
2067 : : contains
2068 : :
2069 : 16429904 : integer function map_symrhg(j1, n1)
2070 : :
2071 : : integer,intent(in) :: j1, n1
2072 : : ! Map into [0,n-1]
2073 : 16429904 : map_symrhg=mod(n1+mod(j1,n1),n1)
2074 : : end function map_symrhg
2075 : :
2076 : : end subroutine symrhg
2077 : : !!***
2078 : :
2079 : : !!****f* m_spacepar/irrzg
2080 : : !! NAME
2081 : : !! irrzg
2082 : : !!
2083 : : !! FUNCTION
2084 : : !! Find the irreducible zone in reciprocal space under the
2085 : : !! symmetry group with real space rotations in symrel(3,3,nsym).
2086 : : !! The (integer) rotation matrices symrel(3,3,nsym) express the new
2087 : : !! real space positions (e.g. rotated atom positions) in REDUCED
2088 : : !! coordinates, i.e. in coordinates expressed as fractions of real space
2089 : : !! primitive translations (atomic coordinates xred). tnons(3,nsym) express
2090 : : !! the associated nonsymmorphic translations, again in reduced coordinates.
2091 : : !! Special data structure created in irrzon.
2092 : : !! First half holds mapping from irr zone to full zone;
2093 : : !! part of second half holds repetition number info.
2094 : : !! work1 is a work array to keep track of grid points found so far.
2095 : : !! In case nspden=2 and nsppol=1, one has to take care of antiferromagnetic
2096 : : !! operations. The subgroup of non-magnetic operations is used
2097 : : !! to generate irrzon(:,:,2) and phnons(:,:,2), while the
2098 : : !! full group is used to generate irrzon(:,:,1) and phnons(:,:,1)
2099 : : !!
2100 : : !! NOTES
2101 : : !! for reference in the near future (2018), some notes: this routine should be duplicated for
2102 : : !! magnetizations in spinorial formalism. The only main difference will
2103 : : !! be that the density is not simply transported to the image point under symrel
2104 : : !! but is a vector which has to be transformed as well
2105 : : !! $S \vec{m} (\vec{x}) = \sigma \vec{m} (S \vec{x} + \tau)$
2106 : : !! $S \vec{m} (\vec{G}) = \sigma exp(+ 2 \pi i \vec{G} \vec{tau}) \vec{m}(S^{-1 T} \vec{G})$
2107 : : !! S is a symop, sigma the AFM sign flip if any, tau the partial non symmorphic translation
2108 : : !! x a position, m the magnetization 3 vector
2109 : : !!
2110 : : !! The phase factor is the same as below for the density, but the collection of elements which
2111 : : !! are equal is more complex: for a 3 or 6 axis the m vector could transform one component
2112 : : !! to a linear combination of several others (I think). Things are not necessarily aligned
2113 : : !! with the z axis.
2114 : : !!
2115 : : !! INPUTS
2116 : : !! nspden=number of spin-density components
2117 : : !! nsppol=1 for unpolarized, 2 for spin-polarized
2118 : : !! nsym=number of symmetry elements in group
2119 : : !! n1,n2,n3=box dimensions of real space grid (or fft box)
2120 : : !! symafm(nsym)=(anti)ferromagnetic part of symmetry operations
2121 : : !! symrel(3,3,nsym)=symmetry matrices in real space (integers)
2122 : : !! tnons(3,nsym)=reduced nonsymmorphic translations
2123 : : !! (symrel and tnons are in terms of real space primitive translations)
2124 : : !!
2125 : : !! OUTPUT
2126 : : !! irrzon(n1*n2*n3,2+(nspden/4),(nspden/nsppol)-3*(nspden/4))=integer array which contains the locations of related
2127 : : !! grid points and the repetition number for each symmetry class.
2128 : : !! phnons(2,n1*n2*n3,(nspden/nsppol)-3*(nspden/4))=phases associated with nonsymmorphic translations
2129 : : !!
2130 : : !! SOURCE
2131 : :
2132 : 6572 : subroutine irrzg(irrzon,nspden,nsppol,nsym,n1,n2,n3,phnons,symafm,symrel,tnons)
2133 : :
2134 : : !Arguments ------------------------------------
2135 : : !scalars
2136 : : integer,intent(in) :: n1,n2,n3,nspden,nsppol,nsym
2137 : : !arrays
2138 : : integer,intent(in) :: symafm(nsym),symrel(3,3,nsym)
2139 : : integer,intent(out) :: irrzon(n1*n2*n3,2,(nspden/nsppol)-3*(nspden/4))
2140 : : real(dp),intent(in) :: tnons(3,nsym)
2141 : : real(dp),intent(out) :: phnons(2,n1*n2*n3,(nspden/nsppol)-3*(nspden/4))
2142 : :
2143 : : !Local variables-------------------------------
2144 : : !scalars
2145 : : integer :: i1,i2,i3,id1,id2,id3,ifft,imagn,ind1,ind2,ipt,irep,isym,izone
2146 : : integer :: izonemax,j1,j2,j3,jj,k1,k2,k3,l1,l2,l3,nfftot,npt,nsym_used
2147 : : integer :: nzone,setzer,sppoldbl
2148 : : real(dp) :: arg,ph1i,ph1r,ph2i,ph2r,tau1,tau2,tau3
2149 : : logical,parameter :: afm_noncoll=.true. ! TRUE if antiferro symmetries are used in non-collinear magnetism
2150 : : character(len=500) :: message
2151 : : !arrays
2152 : 6572 : integer,allocatable :: class(:),iperm(:),symafm_used(:),symrel_used(:,:,:)
2153 : 6572 : integer,allocatable :: work1(:)
2154 : 6572 : real(dp),allocatable :: tnons_used(:,:),work2(:,:)
2155 : : ! *************************************************************************
2156 : :
2157 [ - + + - : 19716 : ABI_MALLOC(class,(nsym))
+ - + - -
+ - + ]
2158 [ - + + - : 13144 : ABI_MALLOC(iperm,(nsym))
- + - + ]
2159 [ - + + - : 19716 : ABI_MALLOC(work1,(n1*n2*n3))
+ - + - -
+ - + ]
2160 [ - + + - : 19716 : ABI_MALLOC(work2,(2,n1*n2*n3))
+ - + - -
+ - + ]
2161 : :
2162 : 6572 : nfftot=n1*n2*n3
2163 : :
2164 : 6572 : id1=n1/2+2
2165 : 6572 : id2=n2/2+2
2166 : 6572 : id3=n3/2+2
2167 : :
2168 [ + + ]: 6572 : sppoldbl=nspden/nsppol;if (nspden==4) sppoldbl=1
2169 : :
2170 [ + + ]: 13175 : do imagn=1,sppoldbl
2171 : :
2172 : : ! Treat in a similar way the case of the full group and the non-magnetic subgroup
2173 : 6603 : nsym_used=0
2174 [ + + ]: 142169 : do isym=1,nsym
2175 [ + + + + : 135566 : if( (imagn==1 .and. sppoldbl==2) .or. symafm(isym)==1 .or. &
+ + ]
2176 : 6603 : & ((nspden==4).and.afm_noncoll) )then
2177 : 134726 : nsym_used=nsym_used+1
2178 : : end if
2179 : : end do
2180 : :
2181 [ + + - + ]: 6603 : if(imagn==2 .and. nsym_used/=nsym/2)then
2182 : : write(message, '(a,a,a,a,a,i4,a,i0)' )&
2183 : 0 : & ' The number of ferromagnetic symmetry operations must be',ch10,&
2184 : 0 : & ' half the total number of operations, while it is observed that',ch10,&
2185 : 0 : & ' nsym=',nsym,' and nsym_magn=',nsym_used
2186 : 0 : ABI_BUG(message)
2187 : : end if
2188 : :
2189 [ - + + - : 19809 : ABI_MALLOC(symafm_used,(nsym_used))
+ - + - -
+ - + -
+ ]
2190 [ - + + - : 19809 : ABI_MALLOC(symrel_used,(3,3,nsym_used))
+ - + - -
+ - + -
+ ]
2191 [ - + + - : 19809 : ABI_MALLOC(tnons_used,(3,nsym_used))
+ - + - -
+ - + -
+ ]
2192 : :
2193 : 6603 : nsym_used=0
2194 [ + + ]: 142169 : do isym=1,nsym
2195 [ + + + + : 135566 : if( (imagn==1 .and. sppoldbl==2) .or. symafm(isym)==1 .or. &
+ + ]
2196 : 6603 : & ((nspden==4).and.afm_noncoll) ) then
2197 : 134726 : nsym_used=nsym_used+1
2198 [ + + + + ]: 1751438 : symrel_used(:,:,nsym_used)=symrel(:,:,isym)
2199 [ + + ]: 538904 : tnons_used(:,nsym_used)=tnons(:,isym)
2200 : 134726 : symafm_used(nsym_used)=symafm(isym)
2201 : : end if
2202 : : end do
2203 [ + + + + ]: 139719 : if ((nspden/=4).or.(.not.afm_noncoll)) symafm_used=1
2204 : :
2205 : :
2206 : : ! Zero out work array--later on, a zero entry will mean that
2207 : : ! a given grid point has not yet been assigned to an ibz point
2208 [ + + ]: 82299287 : work1(1:nfftot)=0
2209 [ + + ]: 82299287 : irrzon(:,2,imagn)=0
2210 : :
2211 : : ! Initialize at G=0 (always in irreducible zone)
2212 : 6603 : nzone=1
2213 : 6603 : irrzon(1,1,imagn)=1
2214 : 6603 : irrzon(1,2,imagn)=nsym_used
2215 : : ! Set phase exp(2*Pi*I*G dot tnons) for G=0
2216 : 6603 : phnons(1,1,imagn)=one
2217 : 6603 : phnons(2,1,imagn)=zero
2218 : 6603 : npt=1
2219 : : ! setting work1(1)=1 indicates that first grid point (G=0) is
2220 : : ! in the iz (irreducible zone)
2221 : 6603 : work1(1)=1
2222 : :
2223 : 6603 : ind1=0
2224 : :
2225 : : ! Loop over reciprocal space grid points:
2226 [ + + ]: 132832 : do i3=1,n3
2227 [ + + ]: 2873533 : do i2=1,n2
2228 [ + + ]: 85159614 : do i1=1,n1
2229 : :
2230 : 82292684 : ind1=ind1+1
2231 : :
2232 : : ! Check to see whether present grid point is equivalent to
2233 : : ! any previously identified ibz point--if not, a new ibz point
2234 : : ! has been found
2235 : :
2236 [ + + ]: 85033385 : if (work1(ind1)==0) then
2237 : :
2238 : : ! A new class has been found.
2239 : :
2240 : : ! Get location of G vector (grid point) centered at 0 0 0
2241 : 10823982 : l3=i3-(i3/id3)*n3-1
2242 : 10823982 : l2=i2-(i2/id2)*n2-1
2243 : 10823982 : l1=i1-(i1/id1)*n1-1
2244 : :
2245 [ + + ]: 117209436 : do isym=1,nsym_used
2246 : :
2247 : : ! Get rotated G vector Gj for each symmetry element
2248 : : ! -- here we use the TRANSPOSE of symrel_used; assuming symrel_used expresses
2249 : : ! the rotation in real space, the transpose is then appropriate
2250 : : ! for G space symmetrization (p. 1172d,e of notes, 2 June 1995).
2251 : : j1=symrel_used(1,1,isym)*l1+&
2252 : 106385454 : & symrel_used(2,1,isym)*l2+symrel_used(3,1,isym)*l3
2253 : : j2=symrel_used(1,2,isym)*l1+&
2254 : 106385454 : & symrel_used(2,2,isym)*l2+symrel_used(3,2,isym)*l3
2255 : : j3=symrel_used(1,3,isym)*l1+&
2256 : 106385454 : & symrel_used(2,3,isym)*l2+symrel_used(3,3,isym)*l3
2257 : :
2258 : : ! Map into [0,n-1]
2259 : 106385454 : k1=mod(n1+mod(j1,n1),n1)
2260 : 106385454 : k2=mod(n2+mod(j2,n2),n2)
2261 : 106385454 : k3=mod(n3+mod(j3,n3),n3)
2262 : :
2263 : : ! Get linear index of rotated point Gj
2264 : 106385454 : ind2=1+k1+n1*(k2+n2*k3)
2265 : :
2266 : : ! Store info for new class:
2267 : 106385454 : class(isym)=ind2
2268 : 106385454 : iperm(isym)=isym
2269 : :
2270 : : ! Setting work array element to 1 indicates grid point has been
2271 : : ! identified with iz point
2272 : 117209436 : work1(ind2)=1
2273 : :
2274 : : ! End of loop on isym
2275 : : end do
2276 : :
2277 : : ! Sort integers into ascending order in each class
2278 : : ! (this lumps together G vectors with the same linear index, i.e.
2279 : : ! groups together symmetries which land on the same Gj)
2280 : 10823982 : call sort_int(nsym_used,class,iperm)
2281 : :
2282 : : ! Check repetition factor (how many identical copies of Gj occur
2283 : : ! from all symmetries applied to G)
2284 : 10823982 : irep=0
2285 [ + + ]: 117209436 : do isym=1,nsym_used
2286 [ + + ]: 117209436 : if (class(isym)==class(1)) then
2287 : 14078455 : irep=irep+1
2288 : : end if
2289 : : end do
2290 : 10823982 : ipt=nsym_used/irep
2291 : :
2292 : : ! Repetition factor must be divisor of nsym_used:
2293 [ - + ]: 10823982 : if (nsym_used/=(ipt*irep)) then
2294 : : write(message, '(a,i5,a,i6,a,a,a,a,a,a)' )&
2295 : 0 : & ' irep=',irep,' not a divisor of nsym_used=',nsym_used,ch10,&
2296 : 0 : & ' This usually indicates that',&
2297 : 0 : & ' the input symmetries do not form a group.',ch10,&
2298 : 0 : & ' Action : check the input symmetries carefully do they',&
2299 : 0 : & ' form a group ? If they do, there is a code bug.'
2300 : 0 : ABI_ERROR(message)
2301 : : end if
2302 : :
2303 : : ! Compute phases for any nonsymmorphic symmetries
2304 : : ! exp(-2*Pi*I*G dot tau(j)) for each symmetry j with
2305 : : ! (possibly zero) nonsymmorphic translation tau(j)
2306 [ + + ]: 117209436 : do jj=1,nsym_used
2307 : : ! First get nonsymmorphic translation and see if nonzero
2308 : : ! (iperm grabs the symmetries in the new order after sorting)
2309 : 106385454 : isym=iperm(jj)
2310 : 106385454 : tau1=tnons_used(1,isym)
2311 : 106385454 : tau2=tnons_used(2,isym)
2312 : 106385454 : tau3=tnons_used(3,isym)
2313 : : if (abs(tau1)>tol12.or.abs(tau2)>tol12&
2314 [ + + + + : 117209436 : & .or.abs(tau3)>tol12) then
+ + ]
2315 : : ! compute exp(-2*Pi*I*G dot tau) using original G
2316 : 17082780 : arg=two_pi*(dble(l1)*tau1+dble(l2)*tau2+dble(l3)*tau3)
2317 : 17082780 : work2(1,jj)=cos(arg)
2318 : 17082780 : work2(2,jj)=-sin(arg)
2319 : : else
2320 : 89302674 : work2(1,jj)=one
2321 : 89302674 : work2(2,jj)=zero
2322 : : end if
2323 : : end do
2324 : :
2325 : : ! All phases arising from symmetries which map to the same
2326 : : ! G vector must actually be the same because
2327 : : ! rho(Strans*G)=exp(2*Pi*I*(G) dot tau_S) rho(G)
2328 : : ! must be satisfied; if exp(2*Pi*I*(G) dot tau_S) can be different
2329 : : ! for two different symmetries S which both take G to the same St*G,
2330 : : ! then the related Fourier components rho(St*G) must VANISH.
2331 : : ! Thus: set "phase" to ZERO here in that case.
2332 : : ! The G mappings occur in sets of irep members; if irep=1 then
2333 : : ! all the G are unique.
2334 : : ! MT 090212:
2335 : : ! In the case of antiferromagn. symetries, one can have
2336 : : ! rho(Strans*G)= -exp(2*Pi*I*(G) dot tau_S) rho(G)
2337 : : ! (look at the minus !)
2338 : : ! A special treatment is then operated on phons.
2339 : : ! The later must be consistent with the use of phnons array
2340 : : ! in symrhg.F90 routine.
2341 : : ! XG 001108 :
2342 : : ! Note that there is a tolerance on the
2343 : : ! accuracy of tnons, especially when they are found from
2344 : : ! the symmetry finder (with xred that might be a bit inaccurate)
2345 [ + + ]: 10823982 : if (irep > 1) then
2346 [ - + + + ]: 18128363 : do jj=1,nsym_used,irep
2347 : 15906979 : setzer=0
2348 : 15906979 : ph1r=work2(1,jj);ph1i=work2(2,jj)
2349 [ + + ]: 55913331 : do j1=jj,jj+irep-1
2350 : 40006352 : ph2r=work2(1,j1);ph2i=work2(2,j1)
2351 [ + + ]: 55913331 : if (((ph2r+ph1r)**2+(ph2i+ph1i)**2) <= tol14) then
2352 [ + + ]: 40026544 : if (setzer/=1) setzer=-1
2353 [ + + ]: 36587394 : else if (((ph2r-ph1r)**2+(ph2i-ph1i)**2) > tol14) then
2354 : 197191 : setzer=1
2355 : : end if
2356 : : end do
2357 : : ! Setzer= 0: phnons are all equal
2358 : : ! Setzer=-1: phnons are equal in absolute value
2359 : : ! Setzer= 1: some phnons are different
2360 [ + + ]: 18128363 : if (setzer/=0) then
2361 [ + + ]: 1675276 : if (setzer==-1) then
2362 [ + + ]: 1574521 : if (afm_noncoll.and.nspden==4) then
2363 : 19484 : arg=symafm_used(iperm(jj))
2364 [ + + + + : 47852 : if (all(symafm_used(iperm(jj:jj+irep-1))==arg)) then
+ + ]
2365 : : setzer=1
2366 : : else
2367 [ + + ]: 41708 : do j1=jj,jj+irep-1
2368 [ + + ]: 98396 : work2(:,j1)=work2(:,j1)*dble(symafm_used(iperm(j1)))
2369 : : end do
2370 : : end if
2371 : : else
2372 : : setzer=1
2373 : : end if
2374 : : end if
2375 [ + + + + : 22943104 : if (setzer==1) work2(:,jj:jj+irep-1)=zero
+ + ]
2376 : : end if
2377 : : end do
2378 : : ! Compress data if irep>1:
2379 : 2221384 : jj=0
2380 [ - + ]: 2221384 : do isym=1,nsym_used,irep
2381 : 15906979 : jj=jj+1
2382 : 15906979 : class(jj)=class(isym)
2383 : 15906979 : work2(1,jj)=work2(1,isym)
2384 [ + + ]: 15906979 : work2(2,jj)=work2(2,isym)
2385 : : end do
2386 : : end if
2387 : :
2388 : : ! Put new unique points into irrzon array:
2389 [ + + ]: 93110063 : irrzon(1+npt:ipt+npt,1,imagn)=class(1:ipt)
2390 : :
2391 : : ! Put repetition number into irrzon array:
2392 : 10823982 : irrzon(1+nzone,2,imagn)=irep
2393 : :
2394 : : ! DEBUG
2395 : : ! write(std_out,'(a,6i7)' )' irrzg : izone,i1,i2,i3,imagn,irrzon(859,2,1)=',&
2396 : : ! & 1+nzone,i1,i2,i3,imagn,irrzon(859,2,1)
2397 : : ! ENDDEBUG
2398 : :
2399 : : ! Put phases (or 0) in phnons array:
2400 [ + + + + ]: 257682225 : phnons(:,1+npt:ipt+npt,imagn)=work2(:,1:ipt)
2401 : :
2402 : : ! Update number of points in irrzon array:
2403 : : ! (irep must divide evenly into nsym_used !)
2404 : 10823982 : npt=npt+ipt
2405 : :
2406 : : ! Update number of classes:
2407 : 10823982 : nzone=nzone+1
2408 : :
2409 : : end if
2410 : : !
2411 : : ! End of loop on reciprocal space points, with indices i1, i2, i3
2412 : : end do
2413 : : end do
2414 : : end do
2415 : :
2416 [ + - ]: 6603 : ABI_SFREE(symafm_used)
2417 [ + - ]: 6603 : ABI_SFREE(symrel_used)
2418 [ + - ]: 13175 : ABI_SFREE(tnons_used)
2419 : :
2420 : : end do ! imagn
2421 : :
2422 : : !Make sure number of real space points accounted for equals actual number of grid points
2423 [ - + ]: 6572 : if (npt/=n1*n2*n3) then
2424 : 0 : write(message, '(a,a,a,a,i10,a,i10,a,a,a,a,a,a,a,a,a)' ) ch10,&
2425 : 0 : & ' irrzg : ERROR -',ch10,&
2426 : 0 : & ' npt=',npt,' and n1*n2*n3=',n1*n2*n3,' are not equal',ch10,&
2427 : 0 : & ' This says that the total of all points in the irreducible',&
2428 : 0 : & ' sector in real space',ch10,&
2429 : 0 : & ' and all symmetrically equivalent',&
2430 : 0 : & ' points, npt, does not equal the actual number',ch10,&
2431 : 0 : & ' of real space grid points.'
2432 : 0 : call wrtout(std_out,message,'COLL')
2433 : : write(message,'(3a)') &
2434 : 0 : & ' This may mean that the input symmetries do not form a group',ch10,&
2435 : 0 : & ' Action : check input symmetries carefully for errors.'
2436 : 0 : ABI_ERROR(message)
2437 : : end if
2438 : :
2439 : : !Perform some checks
2440 [ + + ]: 13175 : do imagn=1,sppoldbl
2441 : :
2442 [ + + ]: 82299287 : do ifft=1,nfftot
2443 [ + - - + ]: 82299287 : if (irrzon(ifft,1,imagn)<1.or.irrzon(ifft,1,imagn)>nfftot) then
2444 : : write(message,'(a,4i0,a,a)')&
2445 : 0 : & ' ifft,irrzon(ifft,1,imagn),nfftot,imagn=',ifft,irrzon(ifft,1,imagn),nfftot,imagn,ch10,&
2446 : 0 : & ' =>irrzon goes outside acceptable bounds.'
2447 : 0 : ABI_BUG(message)
2448 : : end if
2449 : : end do
2450 : :
2451 : 6603 : izonemax=0
2452 [ + + ]: 82305859 : do izone=1,nfftot
2453 : : ! Global bounds
2454 [ + - - + ]: 82292684 : if (irrzon(izone,2,imagn)<0.or.irrzon(izone,2,imagn)>(nsym/imagn)) then
2455 : : write(message, '(a,5i7,a,a)' )&
2456 : 0 : & ' izone,nzone,irrzon(izone,2,imagn),nsym,imagn =',izone,nzone,irrzon(izone,2,imagn),nsym,imagn,ch10,&
2457 : 0 : & ' =>irrzon goes outside acceptable bounds.'
2458 : 0 : ABI_BUG(message)
2459 : : end if
2460 : : ! Second index only goes up to nzone
2461 [ + + ]: 82292684 : if(izonemax==0)then
2462 [ + + ]: 10837188 : if (irrzon(izone,2,imagn)==0)izonemax=izone-1
2463 : : end if
2464 [ + + ]: 10843791 : if(izonemax/=0)then
2465 [ - + ]: 71462099 : if (irrzon(izone,2,imagn)/=0) then
2466 : 0 : message = ' beyond izonemax, irrzon(izone,2,imagn) should be zero'
2467 : 0 : ABI_BUG(message)
2468 : : end if
2469 : : end if
2470 : : end do
2471 : :
2472 : : end do ! imagn
2473 : :
2474 [ - + ]: 6572 : ABI_FREE(class)
2475 [ - + ]: 6572 : ABI_FREE(iperm)
2476 [ - + ]: 6572 : ABI_FREE(work1)
2477 [ - + ]: 6572 : ABI_FREE(work2)
2478 : :
2479 : 6572 : end subroutine irrzg
2480 : : !!***
2481 : :
2482 : : !!****f* m_spacepar/rotate_rho
2483 : : !! NAME
2484 : : !! rotate_rho
2485 : : !!
2486 : : !! FUNCTION
2487 : : !! rotate density in real and reciprocal space
2488 : : !!
2489 : : !! INPUTS
2490 : : !! cplex: if 1, real space 1-order functions on FFT grid are REAL, if 2, COMPLEX
2491 : : !! mpi_enreg=information about MPI parallelization
2492 : : !! nfft=(effective) number of FFT grid points (for this proc) for the "fine" grid (see NOTES in respfn.F90)
2493 : : !! ngfft=array of dimensions for different FFT grids
2494 : : !! nspden=number of spin-density components
2495 : : !! rhor1(cplex*nfft,nspden)=array for Fourier transform of RF electron density
2496 : : !! === if psps%usepaw==1 TODO: extend to PAW
2497 : : !! pawrhoij1(natom) <type(pawrhoij_type)>= 1st-order paw rhoij occupancies and related data
2498 : : !! symrel1=single symmetry operation in real space to apply to rho
2499 : : !! tnon = eventual translation associated to symrel1
2500 : : !!
2501 : : !! OUTPUT
2502 : : !! rhog1_eq(2,nfft)= symmetric density in reciprocal space for equivalent perturbation
2503 : : !! rhor1_eq(cplex*nfft,nspden) = symmetric density in real space for equivalent perturbation
2504 : : !!
2505 : : !! SOURCE
2506 : :
2507 : 25 : subroutine rotate_rho(cplex, itirev, mpi_enreg, nfft, ngfft, nspden, &
2508 : 25 : & rhor1, rhog1_eq, rhor1_eq, symrel1, tnon)
2509 : :
2510 : : !args
2511 : : integer,intent(in) :: cplex, nfft, nspden, itirev
2512 : : integer,intent(in) :: ngfft(18)
2513 : :
2514 : : integer, intent(in) :: symrel1(3,3)
2515 : : real(dp),intent(in) :: tnon(3)
2516 : : real(dp),intent(inout) :: rhor1(cplex*nfft,nspden)
2517 : :
2518 : : real(dp),intent(out) :: rhog1_eq(2,nfft)
2519 : : real(dp),intent(out) :: rhor1_eq(cplex*nfft,nspden)
2520 : :
2521 : : type(MPI_type),intent(in) :: mpi_enreg
2522 : :
2523 : : ! local vars
2524 : : integer :: id1, id2, id3
2525 : : integer :: n1, n2, n3, nd2
2526 : : integer :: l1, l2, l3
2527 : : integer :: i1, i2, i3, ind1, ind2
2528 : : integer :: j1, j2, j3
2529 : : integer :: k1, k2, k3
2530 : : integer :: nproc_fft, ispden, me_fft
2531 : : real(dp) :: arg
2532 : : logical :: t_tnon_nonzero
2533 : :
2534 : : real(dp) :: phnon1(2)
2535 : 25 : real(dp), allocatable :: workg(:,:), workg_eq(:,:)
2536 : : character(len=500) :: message
2537 : : ! *************************************************************************
2538 : :
2539 : 25 : n1=ngfft(1);n2=ngfft(2);n3=ngfft(3);nproc_fft=ngfft(10);me_fft=ngfft(11);nd2=n2/nproc_fft
2540 : :
2541 : 25 : id1=n1/2+2
2542 : 25 : id2=n2/2+2
2543 : 25 : id3=n3/2+2
2544 : :
2545 [ + + + + ]: 86641 : rhog1_eq = zero
2546 [ + + + + ]: 43858 : rhor1_eq = zero
2547 : :
2548 [ - + ]: 25 : if (itirev == 2) then
2549 : 0 : write (message,'(3a,9I4,1a)') 'using time reversal. ',ch10,'Symrel1 = ', symrel1, ch10
2550 : : else
2551 : 25 : write (message,'(3a,9I4,1a)') 'no time reversal. ',ch10,'Symrel1 = ', symrel1, ch10
2552 : : end if
2553 : : !call wrtout(std_out, message, 'COLL')
2554 : :
2555 [ + + + + ]: 82 : t_tnon_nonzero = (any(abs(tnon) > tol12))
2556 : :
2557 : : ! eventually, for FFT parallelization
2558 : : ! call ptabs_fourdp(mpi_enreg,n2,n3,fftn2_distrib,ffti2_local,fftn3_distrib,ffti3_local)
2559 : :
2560 [ - + + - : 75 : ABI_MALLOC(workg,(2,nfft))
+ - + - -
+ - + ]
2561 [ - + + - : 50 : ABI_MALLOC(workg_eq,(2,nfft))
- + - + ]
2562 [ + + ]: 50 : do ispden = 1, nspden
2563 : :
2564 : : ! fft input rhor1 to reciprocal space: uses work* as a buffer
2565 : 25 : call fourdp(cplex,workg,rhor1(:,ispden),-1,mpi_enreg,nfft,1,ngfft,0)
2566 : :
2567 : : ! below taken from irrzg and setsym
2568 : : ! Loop over reciprocal space grid points:
2569 : : ! loop over local points in workg, and get back transform from rhog1,
2570 : : ! which is presumed complete on each proc
2571 : 25 : ind1=0
2572 [ + + ]: 331 : do i3=1,n3
2573 [ + + ]: 3295 : do i2=1,n2
2574 : : ! if(fftn2_distrib(i2)/=me_fft) cycle ! this ind is not to be treated by me_fft
2575 [ + + ]: 32142 : do i1=1,n1
2576 : :
2577 : 28872 : ind1=ind1+1
2578 : : ! r2=ffti2_local(i2+1) - 1
2579 : : ! ind=n1*(nd2*i3+r2)+i1+1 !this is ind in the current proc
2580 : :
2581 : : ! Get location of G vector (grid point) centered at 0 0 0
2582 : 28872 : l1=i1-(i1/id1)*n1-1
2583 : 28872 : l2=i2-(i2/id2)*n2-1
2584 : 28872 : l3=i3-(i3/id3)*n3-1
2585 : :
2586 : : ! Get rotated G vector Gj for each symmetry element
2587 : : ! -- here we use the TRANSPOSE of symrel1; assuming symrel1 expresses
2588 : : ! the rotation in real space, the transpose is then appropriate
2589 : : ! for G space symmetrization (p. 1172d,e of notes, 2 June 1995).
2590 : 28872 : j1=symrel1(1,1)*l1+symrel1(2,1)*l2+symrel1(3,1)*l3
2591 : 28872 : j2=symrel1(1,2)*l1+symrel1(2,2)*l2+symrel1(3,2)*l3
2592 : 28872 : j3=symrel1(1,3)*l1+symrel1(2,3)*l2+symrel1(3,3)*l3
2593 : :
2594 : : ! Map into [0,n-1] and then add 1 for array index in [1,n]
2595 : 28872 : k1=1+mod(n1+mod(j1,n1),n1)
2596 : 28872 : k2=1+mod(n2+mod(j2,n2),n2)
2597 : 28872 : k3=1+mod(n3+mod(j3,n3),n3)
2598 : :
2599 : : ! Get linear index of rotated point Gj
2600 : 28872 : ind2=k1+n1*((k2-1)+n2*(k3-1))
2601 : : ! r2=ffti2_local(j2+1) - 1
2602 : : ! ind=n1*(nd2*j3+r2)+j1+1 !this is ind may be in another proc!!
2603 : :
2604 : 28872 : phnon1(1) = one
2605 : 28872 : phnon1(2) = zero
2606 [ + + ]: 28872 : if (t_tnon_nonzero) then
2607 : : ! compute exp(-2*Pi*I*G dot tau) using original G
2608 : : ! NB: this phase is same as that in irrzg and phnons1, and corresponds to complex conjugate of phase from G to Gj;
2609 : : ! we use it immediately below, to go _to_ workg(ind1)
2610 : : ! TODO : replace this with complex powers of exp(2pi tnon(1)) etc...
2611 : 9600 : arg=two_pi*(dble(l1)*tnon(1)+dble(l2)*tnon(2)+dble(l3)*tnon(3))
2612 : 9600 : phnon1(1) = cos(arg)
2613 : 9600 : phnon1(2) =-sin(arg)
2614 : : end if
2615 : :
2616 : : ! rho(Strans*G)=exp(2*Pi*I*(G) dot tau_S) rho(G)
2617 : : workg_eq (1, ind1) = phnon1(1) * workg(1, ind2) &
2618 : 28872 : & - phnon1(2) * workg(2, ind2)
2619 : : workg_eq (2, ind1) = phnon1(1) * workg(2, ind2) &
2620 : 31836 : & + phnon1(2) * workg(1, ind2)
2621 : :
2622 : : end do
2623 : : end do
2624 : : end do
2625 : :
2626 : : ! accumulate rhog1_eq
2627 [ + - + + : 86641 : if (ispden == 1) rhog1_eq = workg_eq
+ + ]
2628 : :
2629 : : ! FFT back to real space to get rhor1_eq
2630 : : ! Pull out full or spin up density, now symmetrized
2631 : 50 : call fourdp(cplex,workg_eq,rhor1_eq(:,ispden),1,mpi_enreg,nfft,1,ngfft,0)
2632 : :
2633 : : end do !nspden
2634 : :
2635 [ - + ]: 25 : ABI_FREE(workg)
2636 [ - + ]: 25 : ABI_FREE(workg_eq)
2637 : :
2638 : 25 : end subroutine rotate_rho
2639 : : !!***
2640 : :
2641 : : !!****f* m_spacepar/setsym
2642 : : !! NAME
2643 : : !! setsym
2644 : : !!
2645 : : !! FUNCTION
2646 : : !! Set up irreducible zone in G space by direct calculation.
2647 : : !! Do not call this routine if nsym=1 (only identity symmetry).
2648 : : !! Only indsym and symrec get returned if iscf=0.
2649 : : !! symrec needed to symmetrize coordinate gradients in sygrad.
2650 : : !! (symrec is redundant and could be removed later in favor of symrel)
2651 : : !!
2652 : : !! INPUTS
2653 : : !! iscf=(<= 0 =>non-SCF), >0 => SCF
2654 : : !! natom=number of atoms in unit cell
2655 : : !! nfft=(effective) number of FFT grid points (for this processor)
2656 : : !! ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft
2657 : : !! nspden=number of spin-density components
2658 : : !! nsppol=1 for unpolarized, 2 for spin-polarized
2659 : : !! nsym=number of symmetries in space group (at least 1)
2660 : : !! symafm(nsym)=(anti)ferromagnetic part of symmetry operations
2661 : : !! symrel(3,3,nsym)=symmetry operations in terms of real space primitive translations
2662 : : !! tnons(3,nsym)=nonsymmorphic translations of space group in terms
2663 : : !! of real space primitive translations (may be 0)
2664 : : !! typat(natom)=atom type (integer) for each atom
2665 : : !! xred(3,natom)=atomic coordinates in terms of real space primitive translations
2666 : : !!
2667 : : !! OUTPUT
2668 : : !! indsym(4,nsym,natom)=indirect indexing of atom labels--see subroutine
2669 : : !! symatm for definition (if nsym>1)
2670 : : !! irrzon(nfft,2,(nspden/nsppol)-3*(nspden/4))=irreducible zone data
2671 : : !! phnons(2,nfft,(nspden/nsppol)-3*(nspden/4))=nonsymmorphic translation phases
2672 : : !! symrec(3,3,nsym)=symmetry operations in terms of reciprocal
2673 : : !! space primitive translations (if nsym>1)
2674 : : !!
2675 : : !! NOTES
2676 : : !! nsppol and nspden are needed in case of (anti)ferromagnetic symmetry operations
2677 : : !!
2678 : : !! SOURCE
2679 : :
2680 : 9458 : subroutine setsym(indsym,irrzon,iscf,natom,nfft,ngfft,nspden,nsppol,nsym,phnons,&
2681 : 9458 : & symafm,symrec,symrel,tnons,typat,xred)
2682 : :
2683 : : !Arguments ------------------------------------
2684 : : !scalars
2685 : : integer,intent(in) :: iscf,natom,nfft,nspden,nsppol,nsym
2686 : : !arrays
2687 : : integer,intent(in) :: ngfft(18),symafm(nsym),symrel(3,3,nsym),typat(natom)
2688 : : integer,intent(out) :: indsym(4,nsym,natom)
2689 : : integer,intent(inout) :: irrzon(nfft,2,(nspden/nsppol)-3*(nspden/4)) !vz_i
2690 : : integer,intent(out) :: symrec(3,3,nsym)
2691 : : real(dp),intent(in) :: tnons(3,nsym),xred(3,natom)
2692 : : real(dp),intent(out) :: phnons(2,nfft,(nspden/nsppol)-3*(nspden/4))
2693 : :
2694 : : !Local variables-------------------------------
2695 : : !scalars
2696 : : integer :: isym,ierr
2697 : : real(dp) :: tolsym8
2698 : : !arrays
2699 : 9458 : integer,allocatable :: determinant(:)
2700 : : real(dp) :: tsec(2)
2701 : : ! *************************************************************************
2702 : :
2703 : 9458 : call timab(6,1,tsec)
2704 : :
2705 : : !Check that symmetries have unity determinant
2706 [ - + + - : 28374 : ABI_MALLOC(determinant,(nsym))
+ - + - -
+ - + ]
2707 : 9458 : call symdet(determinant,nsym,symrel)
2708 [ - + ]: 9458 : ABI_FREE(determinant)
2709 : :
2710 : :
2711 : : !Get the symmetry matrices in terms of reciprocal basis
2712 [ + + ]: 145978 : do isym=1,nsym
2713 : 145978 : call mati3inv(symrel(:,:,isym),symrec(:,:,isym))
2714 : : end do
2715 : :
2716 : : !Check for group closure
2717 : 9458 : call sg_multable(nsym,symafm,symrel,ierr,tnons=tnons,tnons_tol=tol5)
2718 [ - + ]: 9458 : ABI_CHECK(ierr==0,"Error in group closure")
2719 : :
2720 : 9458 : call sg_multable(nsym,symafm,symrec,ierr)
2721 [ - + ]: 9458 : ABI_CHECK(ierr==0,"Error in group closure")
2722 : :
2723 : : !Obtain a list of rotated atom labels:
2724 : 9458 : tolsym8=tol8
2725 : 9458 : call symatm(indsym,natom,nsym,symrec,tnons,tolsym8,typat,xred,print_indsym=10)
2726 : :
2727 : : !If non-SCF calculation, or nsym==1, do not need IBZ data
2728 [ + + + + ]: 9458 : if ( (iscf>0 .or. iscf==-3) .and. nsym>1 ) then
2729 : : ! Locate irreducible zone in reciprocal space for symmetrization:
2730 : 6106 : call irrzg(irrzon,nspden,nsppol,nsym,ngfft(1),ngfft(2),ngfft(3),phnons,symafm,symrel,tnons)
2731 : : end if
2732 : :
2733 : 9458 : call timab(6,2,tsec)
2734 : :
2735 : : !DEBUG
2736 : : !write(std_out,*)' m_spacepar%setsym : exit '
2737 : : !ENDDEBUG
2738 : :
2739 : 9458 : end subroutine setsym
2740 : : !!***
2741 : :
2742 : : !!****f* ABINIT/hartredq.F90
2743 : : !! NAME
2744 : : !! hartredq.F90
2745 : : !!
2746 : : !! FUNCTION
2747 : : !! Given rho(G), compute the q-gradient of the Hartree potential at q=0
2748 : : !! (=FFT of -rho(G)*G_qdir/pi**2/|G|**4 ) -> Cartesian coordinates
2749 : : !! The calculation is performed in reduced reciprocal space coordinates.
2750 : : !!
2751 : : !! INPUTS
2752 : : !! cplex= if 1, vqgradhartr is REAL, if 2, vqgradhartr is COMPLEX
2753 : : !! gmet(3,3)=metrix tensor in G space in Bohr**-2.
2754 : : !! gprimd(3,3)=reciprocal space dimensional primitive translations
2755 : : !! gsqcut=cutoff value on G**2 for sphere inside fft box.
2756 : : !! (gsqcut=(boxcut**2)*ecut/(2.d0*(Pi**2))
2757 : : !! mpi_enreg=information about MPI parallelization
2758 : : !! nfft=(effective) number of FFT grid points (for this processor)
2759 : : !! ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/input_variables/vargs.htm#ngfft
2760 : : !! qdir= indicates the direction of the q-gradient (1,2 or 3)
2761 : : !! rhog(2,nfft)=electron density in G space
2762 : : !!
2763 : : !! OUTPUT
2764 : : !! vqgradhart(cplex*nfft)=q-gradient of the Hartree potential at q=0in real space, either REAL or COMPLEX
2765 : : !!
2766 : : !! SOURCE
2767 : :
2768 : 591 : subroutine hartredq(cplex,gmet,gsqcut,mpi_enreg,nfft,ngfft,qdir,rhog,vqgradhart)
2769 : :
2770 : : !Arguments ------------------------------------
2771 : : !scalars
2772 : : integer,intent(in) :: cplex,nfft,qdir
2773 : : real(dp),intent(in) :: gsqcut
2774 : : type(MPI_type),intent(in) :: mpi_enreg
2775 : : !arrays
2776 : : integer,intent(in) :: ngfft(18)
2777 : : real(dp),intent(in) :: gmet(3,3),rhog(2,nfft)
2778 : : real(dp),intent(out) :: vqgradhart(cplex*nfft)
2779 : :
2780 : : !Local variables-------------------------------
2781 : : !scalars
2782 : : integer,parameter :: im=2,re=1
2783 : : integer :: i1,i2,i23,i2_local,i3
2784 : : integer :: id1,id2,id3,ig1,ig2,ig3,ii,ii1,me_fft,n1,n2,n3,nproc_fft
2785 : : real(dp) :: cutoff,gfact,gnorm,num
2786 : : real(dp), parameter :: piinv2= piinv*two
2787 : : real(dp),parameter :: tolfix=1.000000001e0_dp
2788 : : !arrays
2789 : 591 : integer, contiguous, pointer :: fftn2_distrib(:),ffti2_local(:)
2790 : 591 : integer, contiguous, pointer :: fftn3_distrib(:),ffti3_local(:)
2791 : 591 : real(dp),allocatable :: work1(:,:)
2792 : : real(dp) :: gvec(3)
2793 : : ! *************************************************************************
2794 : :
2795 : : DBG_ENTER("COLL")
2796 : :
2797 : 591 : n1=ngfft(1); n2=ngfft(2); n3=ngfft(3)
2798 : 591 : nproc_fft = mpi_enreg%nproc_fft; me_fft = mpi_enreg%me_fft
2799 : :
2800 : : !Get the distrib associated with this fft_grid
2801 : 591 : call ptabs_fourdp(mpi_enreg,n2,n3,fftn2_distrib,ffti2_local,fftn3_distrib,ffti3_local)
2802 : :
2803 : : !Initialize a few quantities
2804 : 591 : cutoff=gsqcut*tolfix
2805 [ - + + - : 1773 : ABI_MALLOC(work1,(2,nfft))
+ - + - -
+ - + ]
2806 : 591 : id1=n1/2+2;id2=n2/2+2;id3=n3/2+2
2807 : :
2808 : : !Triple loop on each dimension
2809 [ + + ]: 9894 : do i3=1,n3
2810 : 9303 : ig3=i3-(i3/id3)*n3-1
2811 : :
2812 [ + + ]: 163587 : do i2=1,n2
2813 : 153693 : ig2=i2-(i2/id2)*n2-1
2814 : :
2815 [ + - ]: 162996 : if (fftn2_distrib(i2) == me_fft) then
2816 : 153693 : i2_local = ffti2_local(i2)
2817 : 153693 : i23=n1*(i2_local-1 +(n2/nproc_fft)*(i3-1))
2818 : : !Do the test that eliminates the Gamma point outside of the inner loop
2819 : 153693 : ii1=1
2820 [ + + ]: 153693 : if(i23==0 .and. ig2==0 .and. ig3==0)then
2821 : 591 : ii1=2
2822 : 591 : work1(re,1+i23)=zero
2823 : 591 : work1(im,1+i23)=zero
2824 : : end if
2825 : :
2826 : : ! Final inner loop on the first dimension (note the lower limit)
2827 [ + + ]: 2936625 : do i1=ii1,n1
2828 : 2782932 : ig1=i1-(i1/id1)*n1-1
2829 : 2782932 : ii=i1+i23
2830 : :
2831 [ + + ]: 11131728 : gvec=(/ig1,ig2,ig3/)
2832 : 2782932 : gnorm=normv(gvec,gmet,'r') !'r' is to avoid the 2pi scalling
2833 : :
2834 [ + + ]: 2936625 : if (gnorm**2<=cutoff) then
2835 [ + + ]: 3155120 : num=dot_product(gmet(qdir,:),gvec(:))
2836 : 788780 : gfact=piinv2*num/gnorm**4
2837 : 788780 : work1(re,ii)=-rhog(re,ii)*gfact
2838 : 788780 : work1(im,ii)=-rhog(im,ii)*gfact
2839 : : else
2840 : 1994152 : work1(re,ii)=zero
2841 : 1994152 : work1(im,ii)=zero
2842 : : end if
2843 : :
2844 : : end do ! End loop on i1
2845 : : end if
2846 : :
2847 : : end do ! End loop on i2
2848 : : end do ! End loop on i3
2849 : :
2850 : : ! Fourier Transform the q-gradient of the hartree potential, in reciprocal space it was stored in work1
2851 : 591 : call fourdp(cplex,work1,vqgradhart,1,mpi_enreg,nfft,1,ngfft,0)
2852 : :
2853 [ - + ]: 591 : ABI_FREE(work1)
2854 : :
2855 : : DBG_EXIT("COLL")
2856 : :
2857 : 591 : end subroutine hartredq
2858 : : !!***
2859 : :
2860 : : end module m_spacepar
2861 : : !!***
|
