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 90 : subroutine make_vectornd(cplex,gsqcut,izero,mpi_enreg,natom,nfft,ngfft,nspden,nucdipmom,&
89 90 : & 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 90 : integer,allocatable :: nd_list(:)
113 90 : integer, contiguous, pointer :: fftn2_distrib(:),ffti2_local(:)
114 90 : 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 90 : real(dp),allocatable :: gq(:,:),nd_m(:,:),ndvecr(:),work1(:,:),work2(:,:),work3(:,:)
118 : ! *************************************************************************
119 :
120 90 : 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 90 : nd_atom_tot = 0
126 250 : do iatom = 1, natom
127 464 : if (any(abs(nucdipmom(:,iatom))>tol8)) then
128 96 : 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 90 : call wedge_basis(gprimd,rprimd,rgbasis)
137 :
138 : ! note that nucdipmom is input as vectors in atomic units referenced
139 : ! to cartesian coordinates
140 270 : ABI_MALLOC(nd_list,(nd_atom_tot))
141 270 : ABI_MALLOC(nd_m,(3,nd_atom_tot))
142 250 : nd_atom_tot = 0
143 250 : do iatom = 1, natom
144 464 : if (any(abs(nucdipmom(:,iatom))>tol8)) then
145 96 : nd_atom_tot = nd_atom_tot + 1
146 96 : nd_list(nd_atom_tot) = iatom
147 : ! the following expresses the dipole moment components in units of rprimd translations
148 96 : nd_m(:,nd_atom_tot) = MATMUL(TRANSPOSE(gprimd),nucdipmom(:,iatom))
149 : end if
150 : end do
151 :
152 90 : n1=ngfft(1); n2=ngfft(2); n3=ngfft(3)
153 90 : 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 90 : prefac = -four_pi*j_dpc*FineStructureConstant2/(ucvol*two_pi)
157 :
158 : ! Get the distrib associated with this fft_grid
159 90 : call ptabs_fourdp(mpi_enreg,n2,n3,fftn2_distrib,ffti2_local,fftn3_distrib,ffti3_local)
160 :
161 : ! Initialize a few quantities
162 90 : 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 270 : ABI_MALLOC(gq,(3,max(n1,n2,n3)))
167 360 : do ii=1,3
168 270 : id(ii)=ngfft(ii)/2+2
169 9408 : do ing=1,ngfft(ii)
170 9048 : ig=ing-(ing/id(ii))*ngfft(ii)-1
171 9318 : gq(ii,ing)=ig
172 : end do
173 : end do
174 90 : ig1max=-1;ig2max=-1;ig3max=-1
175 90 : ig1min=n1;ig2min=n2;ig3min=n3
176 :
177 270 : ABI_MALLOC(work1,(2,nfft))
178 180 : ABI_MALLOC(work2,(2,nfft))
179 180 : ABI_MALLOC(work3,(2,nfft))
180 85158234 : work1=zero; work2=zero; work3=zero
181 3106 : id1=n1/2+2;id2=n2/2+2;id3=n3/2+2
182 :
183 : ! Triple loop on each dimension
184 3106 : do i3=1,n3
185 3016 : ig3=i3-(i3/id3)*n3-1
186 : ! Precompute some products that do not depend on i2 and i1
187 3016 : gs3=gq(3,i3)*gq(3,i3)*gmet(3,3)
188 3016 : gqgm23=gq(3,i3)*gmet(2,3)*2
189 3016 : gqgm13=gq(3,i3)*gmet(1,3)*2
190 :
191 146210 : do i2=1,n2
192 143104 : ig2=i2-(i2/id2)*n2-1
193 146120 : if (fftn2_distrib(i2) == me_fft) then
194 143104 : gs2=gs3+ gq(2,i2)*(gq(2,i2)*gmet(2,2)+gqgm23)
195 143104 : gqgm12=gq(2,i2)*gmet(1,2)*2
196 143104 : gqg2p3=gqgm13+gqgm12
197 :
198 143104 : i2_local = ffti2_local(i2)
199 143104 : 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 143104 : 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 9605120 : do i1=ii1,n1
210 9462016 : gs=gs2+ gq(1,i1)*(gq(1,i1)*gmet(1,1)+gqg2p3)
211 9462016 : ig1 = i1 - (i1/id1)*n1 -1
212 9462016 : ii=i1+i23
213 :
214 9462016 : gqred(1) = gq(1,i1); gqred(2) = gq(2,i2); gqred(3) = gq(3,i3)
215 :
216 9605120 : if( (gs .LE. cutoff) .AND. (gs .gt. tol8) )then
217 :
218 9724528 : do iatom = 1, nd_atom_tot
219 5420008 : nd_atom = nd_list(iatom)
220 21680032 : phase = -two_pi*DOT_PRODUCT(xred(:,nd_atom),gqred(:))
221 5420008 : cgr = cmplx(cos(phase),sin(phase))
222 :
223 : ! cross product m x G
224 5420008 : 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 5420008 : mcgc = MATMUL(TRANSPOSE(gprimd),mcgc)
231 :
232 5420008 : work1(re,ii) = work1(re,ii) + real(prefac*cgr*mcgc(1)/gs)
233 5420008 : work2(re,ii) = work2(re,ii) + real(prefac*cgr*mcgc(2)/gs)
234 5420008 : work3(re,ii) = work3(re,ii) + real(prefac*cgr*mcgc(3)/gs)
235 :
236 5420008 : work1(im,ii) = work1(im,ii) + aimag(prefac*cgr*mcgc(1)/gs)
237 5420008 : work2(im,ii) = work2(im,ii) + aimag(prefac*cgr*mcgc(2)/gs)
238 9724528 : work3(im,ii) = work3(im,ii) + aimag(prefac*cgr*mcgc(3)/gs)
239 :
240 : end do
241 : else
242 : ! gs>cutoff
243 5157496 : work1(re,ii)=zero
244 5157496 : work1(im,ii)=zero
245 5157496 : work2(re,ii)=zero
246 5157496 : work2(im,ii)=zero
247 5157496 : work3(re,ii)=zero
248 5157496 : 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 90 : ABI_FREE(gq)
257 90 : ABI_FREE(nd_list)
258 90 : ABI_FREE(nd_m)
259 :
260 90 : if ( izero .EQ. 1 ) then
261 : ! Set contribution of unbalanced components to zero
262 :
263 85 : call zerosym(work1,2,n1,n2,n3,comm_fft=mpi_enreg%comm_fft,distribfft=mpi_enreg%distribfft)
264 85 : call zerosym(work2,2,n1,n2,n3,comm_fft=mpi_enreg%comm_fft,distribfft=mpi_enreg%distribfft)
265 85 : 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 28635528 : vectornd = zero
272 : ! Fourier Transform
273 270 : ABI_MALLOC(ndvecr,(cplex*nfft))
274 9462106 : ndvecr=zero
275 90 : call fourdp(cplex,work1,ndvecr,1,mpi_enreg,nfft,1,ngfft,0)
276 9462106 : vectornd(:,1,1)=ndvecr(:)
277 83034 : if (nspden .GE. 2) vectornd(:,2,1) = ndvecr(:)
278 90 : ABI_FREE(work1)
279 :
280 9462106 : ndvecr=zero
281 90 : call fourdp(cplex,work2,ndvecr,1,mpi_enreg,nfft,1,ngfft,0)
282 9462106 : vectornd(:,1,2) = ndvecr(:)
283 83034 : if (nspden .GE. 2) vectornd(:,2,2) = ndvecr(:)
284 90 : ABI_FREE(work2)
285 :
286 9462106 : ndvecr=zero
287 90 : call fourdp(cplex,work3,ndvecr,1,mpi_enreg,nfft,1,ngfft,0)
288 9462106 : vectornd(:,1,3) = ndvecr(:)
289 83034 : if (nspden .GE. 2) vectornd(:,2,3) = ndvecr(:)
290 90 : ABI_FREE(work3)
291 90 : ABI_FREE(ndvecr)
292 :
293 90 : 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 83430 : subroutine hartre(cplex,gsqcut,icutcoul,izero,mpi_enreg,nfft,ngfft,nkpt,&
551 83430 : 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 83430 : integer, contiguous, pointer :: fftn2_distrib(:),ffti2_local(:)
578 83430 : integer, contiguous, pointer :: fftn3_distrib(:),ffti3_local(:)
579 : real(dp) :: gmet(3,3),gprimd(3,3),qpt_(3),rmet(3,3),tsec(2)
580 83430 : real(dp),allocatable :: gcutoff(:)
581 83430 : real(dp),allocatable :: gq(:,:),work1(:,:)
582 : ! *************************************************************************
583 :
584 : ! Keep track of total time spent in hartre
585 83430 : call timab(10,1,tsec)
586 :
587 : ! Check that cplex has an allowed value
588 83430 : 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 83430 : call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
595 :
596 83430 : n1=ngfft(1); n2=ngfft(2); n3=ngfft(3)
597 83430 : nproc_fft = mpi_enreg%nproc_fft; me_fft = mpi_enreg%me_fft
598 :
599 : ! Get the distrib associated with this fft_grid
600 83430 : call ptabs_fourdp(mpi_enreg,n2,n3,fftn2_distrib,ffti2_local,fftn3_distrib,ffti3_local)
601 :
602 : ! Initialize a few quantities
603 83430 : cutoff=gsqcut*tolfix
604 83430 : if(present(qpt))then
605 21014 : qpt_=qpt
606 : else
607 62416 : qpt_=zero
608 : end if
609 83430 : qeq0=0
610 83430 : 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 83430 : 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 83430 : 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 83430 : nog0=0; if (qeq0==1 .or. icutcoul==55) nog0=1
634 :
635 : !Initialize Gcut-off array from m_gtermcutoff
636 83430 : 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 250290 : ABI_MALLOC(gq,(3,max(n1,n2,n3)))
641 333720 : do ii=1,3
642 250290 : id(ii)=ngfft(ii)/2+2
643 5955738 : do ing=1,ngfft(ii)
644 5622018 : ig=ing-(ing/id(ii))*ngfft(ii)-1
645 5872308 : gq(ii,ing)=ig+qpt_(ii)
646 : end do
647 : end do
648 83430 : ig1max=-1;ig2max=-1;ig3max=-1
649 83430 : ig1min=n1;ig2min=n2;ig3min=n3
650 :
651 250290 : ABI_MALLOC(work1,(2,nfft))
652 83430 : id1=n1/2+2;id2=n2/2+2;id3=n3/2+2
653 :
654 : ! Triple loop on each dimension
655 2012004 : do i3=1,n3
656 1928574 : ig3=i3-(i3/id3)*n3-1
657 : ! Precompute some products that do not depend on i2 and i1
658 1928574 : gs3=gq(3,i3)*gq(3,i3)*gmet(3,3)
659 1928574 : gqgm23=gq(3,i3)*gmet(2,3)*2
660 1928574 : gqgm13=gq(3,i3)*gmet(1,3)*2
661 :
662 52982892 : do i2=1,n2
663 50970888 : ig2=i2-(i2/id2)*n2-1
664 52899462 : if (fftn2_distrib(i2) == me_fft) then
665 49139912 : gs2=gs3+ gq(2,i2)*(gq(2,i2)*gmet(2,2)+gqgm23)
666 49139912 : gqgm12=gq(2,i2)*gmet(1,2)*2
667 49139912 : gqg2p3=gqgm13+gqgm12
668 :
669 49139912 : i2_local = ffti2_local(i2)
670 49139912 : 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 49139912 : ii1=1
673 : ! if(i23==0 .and. qeq0==1 .and. ig2==0 .and. ig3==0)then
674 49139912 : if(i23==0 .and. nog0==1 .and. ig2==0 .and. ig3==0)then
675 75624 : ii1=2
676 75624 : work1(re,1+i23)=zero
677 75624 : work1(im,1+i23)=zero
678 : end if
679 :
680 : ! Final inner loop on the first dimension (note the lower limit)
681 1939638934 : do i1=ii1,n1
682 1890499022 : gs=gs2+ gq(1,i1)*(gq(1,i1)*gmet(1,1)+gqg2p3)
683 1890499022 : ii=i1+i23
684 :
685 1939638934 : if(gs<=cutoff)then
686 : ! Identify min/max indexes (to cancel unbalanced contributions later)
687 : ! Count (q+g)-vectors with similar norm
688 732578038 : 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 732578038 : den=piinv/gs*gcutoff(ii)
696 732578038 : work1(re,ii)=rhog(re,ii)*den
697 732578038 : work1(im,ii)=rhog(im,ii)*den
698 : else
699 : ! gs>cutoff
700 1157920984 : work1(re,ii)=zero
701 1157920984 : 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 83430 : ABI_FREE(gq)
710 :
711 83430 : if (izero==1) then
712 : ! Set contribution of unbalanced components to zero
713 :
714 22858 : if (qeq0==1) then !q=0
715 22858 : 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 83430 : call fourdp(cplex,work1,vhartr,1,mpi_enreg,nfft,1,ngfft,0)
741 :
742 83430 : ABI_FREE(gcutoff)
743 83430 : ABI_FREE(work1)
744 :
745 83430 : call timab(10,2,tsec)
746 :
747 83430 : 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 5081469 : 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 5081469 : 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 5081469 : 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 5081469 : me_g0 = mpi_enreg%me_g0
811 :
812 5081469 : ar=zero
813 5081469 : if(present(ar_im)) ar_im=zero
814 :
815 5081469 : l_gpu_thread_limit=0; if(present(gpu_thread_limit)) l_gpu_thread_limit=gpu_thread_limit
816 :
817 3927669 : 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 5081469 : if(istwf_k==1)then
824 :
825 : ! No filter
826 4912887 : if(filter==0)then
827 : !$OMP PARALLEL DO REDUCTION(+:ar)
828 988645363 : do ipw=1,npw
829 988645363 : ar=ar+diag(ipw)*(vect(1,ipw)*vect1(1,ipw)+vect(2,ipw)*vect1(2,ipw))
830 : end do
831 4861559 : if(nspinor==2)then
832 : !$OMP PARALLEL DO REDUCTION(+:ar) PRIVATE(jpw)
833 53472556 : do ipw=1+npw,2*npw
834 53194216 : jpw=ipw-npw
835 53472556 : ar=ar+diag(jpw)*(vect(1,ipw)*vect1(1,ipw)+vect(2,ipw)*vect1(2,ipw))
836 : end do
837 : end if
838 4861559 : if(use_ndo==1)then
839 : !$OMP PARALLEL DO REDUCTION(+:ar_im)
840 191520106 : do ipw=1,npw
841 191520106 : ar_im=ar_im+diag(ipw)*(vect1(1,ipw)*vect(2,ipw)-vect1(2,ipw)*vect(1,ipw))
842 : end do
843 510096 : if(nspinor == 2) then
844 : !$OMP PARALLEL DO REDUCTION(+:ar_im) PRIVATE(jpw)
845 714420 : do ipw=1+npw,2*npw
846 710892 : jpw=ipw-npw
847 714420 : 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 11366416 : do ipw=1,npw
857 11366416 : if(diag(ipw)<hugevalue)then
858 10966672 : ar=ar+diag(ipw)*(vect(1,ipw)*vect1(1,ipw)+vect(2,ipw)*vect1(2,ipw))
859 : end if
860 : end do
861 51328 : 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 51328 : 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 168582 : else if(istwf_k>=2)then
894 :
895 168582 : if(filter==0)then
896 168582 : i1=1
897 168582 : if(istwf_k==2 .and. me_g0==1)then ! MPIWF need to know which proc has G=0
898 72567 : ar=half*diag(1)*vect(1,1)*vect1(1,1) ; i1=2
899 : end if
900 :
901 : !$OMP PARALLEL DO REDUCTION(+:ar)
902 52850145 : do ipw=i1,npw
903 52850145 : 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 168582 : ar=two*ar
924 :
925 : end if ! istwf_k
926 :
927 : !MPIWF need to make reduction on ar and ai .
928 5081469 : 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 5081469 : if (l_gpu_thread_limit /= 0) call xomp_set_num_threads(nthreads_bak)
934 :
935 5081469 : 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 2124877 : subroutine meanvalue_g_batch(ar, diag, filter, istwf_k, mpi_enreg, npw, nspinor, ndat, &
969 2124877 : 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 2124877 : 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 2124877 : l_gpu_option = ABI_GPU_DISABLED; if (present(gpu_option)) l_gpu_option = gpu_option
998 2124877 : l_gpu_thread_limit=0; if(present(gpu_thread_limit)) l_gpu_thread_limit=gpu_thread_limit
999 2124877 : if(l_gpu_option==ABI_GPU_OPENMP) l_gpu_thread_limit=0
1000 :
1001 4712690 : ar(:) = zero
1002 :
1003 2124877 : 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 4124997 : do idat = 1, ndat
1010 : local_ar = zero
1011 : !$OMP PARALLEL DO REDUCTION(+:local_ar) PRIVATE(jpw)
1012 648016004 : do ipw = 1, npw*nspinor
1013 645818316 : jpw = mod(ipw-1, npw) + 1
1014 : local_ar = local_ar + diag(jpw) * ( vect(1, ipw+(idat-1)*npw*nspinor)**2 &
1015 648016004 : +vect(2, ipw+(idat-1)*npw*nspinor)**2)
1016 : end do
1017 4124997 : ar(idat) = local_ar
1018 : end do
1019 :
1020 1927309 : if (mpi_enreg%paral_kgb == 1) then
1021 110628 : call xmpi_sum(ar, ndat, mpi_enreg%comm_bandspinorfft, ierr)
1022 : end if
1023 :
1024 197568 : else if (istwf_k >= 2 .and. filter == 0 .and. use_ndo == 0) then
1025 :
1026 197568 : me_g0 = mpi_enreg%me_g0
1027 197568 : i1 = 1
1028 197568 : 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 587693 : do idat = 1, ndat
1035 390125 : local_ar = zero
1036 390125 : if (i1 == 2) then
1037 111033 : local_ar = half * diag(1) * vect(1, 1+(idat-1)*npw)**2
1038 : end if
1039 : !$OMP PARALLEL DO REDUCTION(+:local_ar)
1040 146393962 : do ipw = i1, npw
1041 : local_ar = local_ar + diag(ipw) * (vect(1, ipw+(idat-1)*npw)**2 &
1042 146393962 : +vect(2, ipw+(idat-1)*npw)**2)
1043 : end do
1044 587693 : ar(idat) = two * local_ar
1045 : end do
1046 :
1047 197568 : 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 2124877 : if (l_gpu_thread_limit /= 0) call xomp_set_num_threads(nthreads_bak)
1076 :
1077 2124877 : 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 71085 : subroutine symrhg(cplex,gprimd,irrzon,mpi_enreg,nfft,nfftot,ngfft,nspden,nsppol,nsym,&
1612 71085 : & 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 71085 : integer,allocatable :: isymg(:)
1638 71085 : integer, contiguous, pointer :: fftn2_distrib(:),ffti2_local(:)
1639 71085 : integer, contiguous, pointer :: fftn3_distrib(:),ffti3_local(:)
1640 : real(dp) :: tsec(2)
1641 71085 : real(dp),allocatable :: magngx(:,:),magngy(:,:),magngz(:,:)
1642 71085 : real(dp),allocatable :: rhosu1_arr(:),rhosu2_arr(:),work(:)
1643 71085 : 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 213255 : ABI_MALLOC(work,(cplex*nfft))
1649 :
1650 : !Special treatment for spin-polarized case
1651 71085 : 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 6505 : call timab(17,1,tsec)
1655 112434012 : work(:)=rhor(:,1) ! up => work
1656 112434012 : rhor(:,1)=rhor(:,1)+rhor(:,2) ! up+down
1657 112434012 : rhor(:,2)=work(:) ! work => up
1658 6505 : call timab(17,2,tsec)
1659 : end if
1660 :
1661 : !Special treatment for antiferromagnetism case
1662 71085 : 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 71085 : if(nspden==4) then
1674 2809 : call timab(17,1,tsec)
1675 : !FR the half factors missing are recovered in dfpt_mkvxc_noncoll and dfpt_accrho
1676 24733921 : rhor(:,1)=rhor(:,1)+rhor(:,4) !nup+ndown
1677 24733921 : rhor(:,2)=rhor(:,2)-rhor(:,1) !mx (n+mx-n)
1678 24733921 : rhor(:,3)=rhor(:,3)-rhor(:,1) !my (n+my-n)
1679 24733921 : rhor(:,4)=rhor(:,1)-two*rhor(:,4) !mz=n-2ndown
1680 2809 : call timab(17,2,tsec)
1681 : end if
1682 :
1683 :
1684 71085 : if(nsym==1)then
1685 :
1686 25831 : 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 347233435 : work(:)=rhor(:,1)
1693 25831 : 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 45254 : nspden_eff=nspden;if (nspden==4) nspden_eff=1
1700 95253 : do ispden=nspden_eff,1,-1
1701 :
1702 : ! Prepare the density to be symmetrized, in the reciprocal space
1703 49999 : if(nspden==1 .or. nsppol==2 .or. (nspden==4.and.(.not.afm_noncoll)))then
1704 1083779 : imagn=1
1705 : nsym_used=0
1706 1083779 : do isym=1,nsym
1707 1083779 : 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 2660 : else if(nspden==2 .and. nsppol==1)then ! antiferromagnetic case
1713 846 : imagn=ispden
1714 846 : nsym_used=nsym/ispden
1715 1814 : else if (nspden==4) then
1716 1814 : imagn=1
1717 1814 : 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 745284277 : work(:)=rhor(:,ispden)
1725 49999 : call fourdp(cplex,rhog,work,-1,mpi_enreg,nfft,1,ngfft,0)
1726 49999 : if (nspden==4) then
1727 5442 : ABI_MALLOC(magngx,(2,nfft))
1728 3628 : ABI_MALLOC(magngy,(2,nfft))
1729 3628 : ABI_MALLOC(magngz,(2,nfft))
1730 13723743 : work(:)=rhor(:,2)
1731 1814 : call fourdp(cplex,magngx,work,-1,mpi_enreg,nfft,1,ngfft,0)
1732 13723743 : work(:)=rhor(:,3)
1733 1814 : call fourdp(cplex,magngy,work,-1,mpi_enreg,nfft,1,ngfft,0)
1734 13723743 : work(:)=rhor(:,4)
1735 1814 : 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 49999 : call timab(17,1,tsec)
1740 :
1741 49999 : 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 49999 : 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 750733068 : izone_max=count(irrzon(:,2,imagn)>0)
1751 149997 : ABI_MALLOC(rhosu1_arr,(izone_max))
1752 99998 : ABI_MALLOC(rhosu2_arr,(izone_max))
1753 :
1754 49999 : numpt=0
1755 109301315 : do izone=1,nfftot
1756 :
1757 : ! Get repetition number
1758 109301315 : rep=irrzon(izone,2,imagn)
1759 109301315 : if(rep==0)exit
1760 :
1761 : ! Compute number of unique points in this symm class:
1762 109251316 : nup=nsym_used/rep
1763 :
1764 : ! Accumulate charge over equivalent points
1765 109251316 : rhosu1=zero
1766 109251316 : rhosu2=zero
1767 859754234 : do iup=1,nup
1768 750502918 : ind=irrzon(iup+numpt,1,imagn)
1769 750502918 : j=ind-1;j1=modulo(j,n1);j2=modulo(j/n1,n2);j3=j/(n1*n2);
1770 859754234 : if(fftn2_distrib(j2+1)==me_fft) then ! this ind is to be treated by me_fft
1771 734573494 : r2=ffti2_local(j2+1) - 1
1772 734573494 : 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 734573494 : & -rhog(2,ind)*phnons(2,iup+numpt,imagn)
1775 : rhosu2=rhosu2+rhog(2,ind)*phnons(1,iup+numpt,imagn)&
1776 734573494 : & +rhog(1,ind)*phnons(2,iup+numpt,imagn)
1777 : end if
1778 :
1779 : end do
1780 109251316 : rhosu1=rhosu1/dble(nup)
1781 109251316 : rhosu2=rhosu2/dble(nup)
1782 109251316 : rhosu1_arr(izone)=rhosu1
1783 109251316 : rhosu2_arr(izone)=rhosu2
1784 : ! Keep index of how many points have been considered:
1785 109301315 : numpt=numpt+nup
1786 :
1787 : end do ! End loop over izone
1788 :
1789 : ! Reduction in case of FFT parallelization
1790 49999 : if(mpi_enreg%nproc_fft>1)then
1791 2640 : spaceComm=mpi_enreg%comm_fft
1792 2640 : call xmpi_sum(rhosu1_arr,spaceComm,ier)
1793 2640 : call xmpi_sum(rhosu2_arr,spaceComm,ier)
1794 : end if
1795 :
1796 : ! Now symmetrize the density
1797 : numpt=0
1798 109301315 : do izone=1,nfftot
1799 :
1800 : ! Get repetition number
1801 109301315 : rep=irrzon(izone,2,imagn)
1802 109301315 : if(rep==0)exit
1803 :
1804 : ! Compute number of unique points in this symm class:
1805 109251316 : nup=nsym_used/rep
1806 :
1807 : ! Define symmetrized rho(G) at equivalent points:
1808 859754234 : do iup=1,nup
1809 750502918 : ind=irrzon(iup+numpt,1,imagn)
1810 : ! decompose ind-1=n1(n2 j3+ j2)+j1
1811 750502918 : j=ind-1;j1=modulo(j,n1);j2=modulo(j/n1,n2);j3=j/(n1*n2);
1812 859754234 : if(fftn2_distrib(j2+1)==me_fft) then ! this ind is to be treated by me_fft
1813 734573494 : r2=ffti2_local(j2+1) - 1
1814 : ! ind in the proc ind-1=n1(nd2 j3+ r2)+j1
1815 734573494 : 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 734573494 : & +rhosu2_arr(izone)*phnons(2,iup+numpt,imagn)
1818 : rhog(2,ind)=rhosu2_arr(izone)*phnons(1,iup+numpt,imagn)&
1819 734573494 : & -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 109301315 : numpt=numpt+nup
1825 :
1826 : end do ! End loop over izone
1827 :
1828 49999 : ABI_FREE(rhosu1_arr)
1829 49999 : ABI_FREE(rhosu2_arr)
1830 :
1831 : ! The following is only valid for magnetization
1832 : ! ---------------------------------------------
1833 49999 : if (nspden==4) then
1834 :
1835 1814 : id1=n1/2+2
1836 1814 : id2=n2/2+2
1837 1814 : id3=n3/2+2
1838 :
1839 : ! Transfer symmetries in cartesian coordinates
1840 : ! Compute symmetries in reciprocal space in cartesian coordinates
1841 5442 : ABI_MALLOC(symrec_cart,(3,3,nsym_used))
1842 3628 : ABI_MALLOC(symrel_cart,(3,3,nsym_used))
1843 5442 : ABI_MALLOC(symafm_used,(nsym_used))
1844 5442 : ABI_MALLOC(tnons_used,(3,nsym_used))
1845 3628 : ABI_MALLOC(sym_det,(nsym_used))
1846 13668 : jsym=0
1847 13668 : do isym=1,nsym
1848 : if (symafm(isym)/=1.and.(.not.afm_noncoll)) cycle
1849 11854 : jsym=jsym+1
1850 47416 : tnons_used(:,jsym)=tnons(:,isym)
1851 11854 : symafm_used(jsym)=dble(symafm(isym))
1852 11854 : call symredcart(rprimd,gprimd,symrel_cart(:,:,jsym),symrel(:,:,isym))
1853 11854 : 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 13668 : & symrel_cart(3,2,isym)*symrel_cart(2,3,isym)*symrel_cart(1,1,isym))
1860 : end do
1861 :
1862 13830239 : numpt=count(irrzon(:,1,imagn)>0)
1863 5442 : ABI_MALLOC(isymg,(numpt))
1864 13830239 : isymg=0
1865 5442 : ABI_MALLOC(rhosu1_arr,(3*izone_max))
1866 3628 : 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 1814 : numpt=0
1872 5097073 : do izone=1,izone_max
1873 5095259 : magxsu1=zero;magxsu2=zero
1874 5095259 : magysu1=zero;magysu2=zero
1875 5095259 : magzsu1=zero;magzsu2=zero
1876 5095259 : ind=irrzon(1+numpt,1,1)
1877 5095259 : rep=irrzon(izone,2,1)
1878 5095259 : nup=nsym_used/rep
1879 : ! Get coordinates in the range [0,n-1]
1880 5095259 : 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 5095259 : l3=l3-(l3/id3)*n3
1884 5095259 : l2=l2-(l2/id2)*n2
1885 5095259 : l1=l1-(l1/id1)*n1
1886 :
1887 5095259 : jsym=0
1888 20453115 : do isym=1,nsym
1889 : if (symafm(isym)/=1.and.(.not.afm_noncoll)) cycle
1890 15357856 : 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 15357856 : j1=symrel(1,1,isym)*l1+symrel(2,1,isym)*l2+symrel(3,1,isym)*l3
1894 15357856 : j2=symrel(1,2,isym)*l1+symrel(2,2,isym)*l2+symrel(3,2,isym)*l3
1895 15357856 : j3=symrel(1,3,isym)*l1+symrel(2,3,isym)*l2+symrel(3,3,isym)*l3
1896 15357856 : k1=map_symrhg(j1,n1);k2=map_symrhg(j2,n2);k3=map_symrhg(j3,n3)
1897 15357856 : indsy=1+k1+n1*(k2+n2*k3)
1898 15357856 : ind2=-1;iup=numpt
1899 57490340 : do while (ind2/=indsy.and.iup<numpt+nup)
1900 42132484 : iup=iup+1;ind2=irrzon(iup,1,1)
1901 : end do
1902 15357856 : if (ind2/=indsy) then
1903 0 : ABI_ERROR("ind2/=indsy in symrhg !")
1904 : end if
1905 15357856 : if (isymg(iup)==0) isymg(iup)=jsym
1906 20453115 : if(fftn2_distrib(modulo((indsy-1)/n1,n2) + 1) == me_fft ) then ! this is indsy is to be treated by me_fft
1907 15245536 : 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 15245536 : tau1=tnons_used(1,jsym)
1911 15245536 : tau2=tnons_used(2,jsym)
1912 15245536 : tau3=tnons_used(3,jsym)
1913 15245536 : 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 1249408 : arg=two_pi*(dble(l1)*tau1+dble(l2)*tau2+dble(l3)*tau3)
1916 1249408 : phr=cos(arg)
1917 1249408 : phi=-sin(arg)
1918 : else
1919 : phr=one
1920 : phi=zero
1921 : end if
1922 15245536 : phr=phr*symafm_used(jsym)
1923 15245536 : 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 15245536 : 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 15245536 : 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 15245536 : 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 15245536 : 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 15245536 : 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 15245536 : 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 15245536 : magxsu1=magxsu1+mxr*phr-mxi*phi;magxsu2=magxsu2+mxi*phr+mxr*phi
1946 15245536 : magysu1=magysu1+myr*phr-myi*phi;magysu2=magysu2+myi*phr+myr*phi
1947 15245536 : magzsu1=magzsu1+mzr*phr-mzi*phi;magzsu2=magzsu2+mzi*phr+mzr*phi
1948 : end if
1949 : end do
1950 5095259 : rhosu1_arr(3*izone-2)=magxsu1/dble(nsym_used)
1951 5095259 : rhosu1_arr(3*izone-1)=magysu1/dble(nsym_used)
1952 5095259 : rhosu1_arr(3*izone )=magzsu1/dble(nsym_used)
1953 5095259 : rhosu2_arr(3*izone-2)=magxsu2/dble(nsym_used)
1954 5095259 : rhosu2_arr(3*izone-1)=magysu2/dble(nsym_used)
1955 5095259 : rhosu2_arr(3*izone )=magzsu2/dble(nsym_used)
1956 5097073 : numpt=numpt+nup
1957 : end do
1958 :
1959 : ! Reduction in case of FFT parallelization
1960 1814 : if(mpi_enreg%nproc_fft>1)then
1961 52 : spaceComm=mpi_enreg%comm_fft
1962 52 : call xmpi_sum(rhosu1_arr,spaceComm,ier)
1963 52 : 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 5097073 : do izone=1,izone_max
1970 5095259 : rep=irrzon(izone,2,imagn)
1971 5095259 : nup=nsym_used/rep
1972 18923684 : do iup=1,nup
1973 13828425 : ind=irrzon(iup+numpt,1,imagn)
1974 : ! Get coordinates in the range [0,n-1]
1975 13828425 : 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 13828425 : l3=j3-(j3/id3)*n3
1979 13828425 : l2=j2-(j2/id2)*n2
1980 13828425 : l1=j1-(j1/id1)*n1
1981 18923684 : if(fftn2_distrib(j2+1)==me_fft) then ! this ind is to be treated by me_fft
1982 13721929 : r2=ffti2_local(j2+1) - 1
1983 13721929 : ind=n1*(nd2*j3+r2)+j1+1 ! this is ind in the current proc
1984 13721929 : jsym=isymg(iup+numpt)
1985 13721929 : if (jsym==0) then
1986 0 : ABI_ERROR("jsym=0 in symrhg !")
1987 : end if
1988 13721929 : magxsu1=rhosu1_arr(3*izone-2);magxsu2=rhosu2_arr(3*izone-2)
1989 13721929 : magysu1=rhosu1_arr(3*izone-1);magysu2=rhosu2_arr(3*izone-1)
1990 13721929 : 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 13721929 : tau1=tnons_used(1,jsym)
1993 13721929 : tau2=tnons_used(2,jsym)
1994 13721929 : tau3=tnons_used(3,jsym)
1995 13721929 : 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 848614 : arg=two_pi*(dble(l1)*tau1+dble(l2)*tau2+dble(l3)*tau3)
1998 848614 : phr=cos(arg)
1999 848614 : phi=-sin(arg)
2000 : else
2001 : phr=one
2002 : phi=zero
2003 : end if
2004 13721929 : phr=phr*symafm_used(jsym)
2005 13721929 : 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 13721929 : mxr=sym_det(jsym)*(symrec_cart(1,1,jsym)*magxsu1+symrec_cart(2,1,jsym)*magysu1+symrec_cart(3,1,jsym)*magzsu1)
2012 13721929 : mxi=sym_det(jsym)*(symrec_cart(1,1,jsym)*magxsu2+symrec_cart(2,1,jsym)*magysu2+symrec_cart(3,1,jsym)*magzsu2)
2013 13721929 : myr=sym_det(jsym)*(symrec_cart(1,2,jsym)*magxsu1+symrec_cart(2,2,jsym)*magysu1+symrec_cart(3,2,jsym)*magzsu1)
2014 13721929 : myi=sym_det(jsym)*(symrec_cart(1,2,jsym)*magxsu2+symrec_cart(2,2,jsym)*magysu2+symrec_cart(3,2,jsym)*magzsu2)
2015 13721929 : mzr=sym_det(jsym)*(symrec_cart(1,3,jsym)*magxsu1+symrec_cart(2,3,jsym)*magysu1+symrec_cart(3,3,jsym)*magzsu1)
2016 13721929 : 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 13721929 : magngx(1,ind)=mxr*phr-mxi*phi
2024 13721929 : magngx(2,ind)=mxi*phr+mxr*phi
2025 13721929 : magngy(1,ind)=myr*phr-myi*phi
2026 13721929 : magngy(2,ind)=myi*phr+myr*phi
2027 13721929 : magngz(1,ind)=mzr*phr-mzi*phi
2028 13721929 : magngz(2,ind)=mzi*phr+mzr*phi
2029 : end if
2030 : end do
2031 5097073 : numpt=numpt+nup
2032 : end do
2033 1814 : ABI_FREE(isymg)
2034 1814 : ABI_FREE(rhosu1_arr)
2035 1814 : ABI_FREE(rhosu2_arr)
2036 1814 : ABI_FREE(symrec_cart)
2037 1814 : ABI_FREE(symrel_cart)
2038 1814 : ABI_FREE(sym_det)
2039 1814 : ABI_FREE(symafm_used)
2040 1814 : ABI_FREE(tnons_used)
2041 :
2042 : end if ! nspden==4
2043 :
2044 49999 : call timab(17,2,tsec)
2045 :
2046 : ! Pull out full or spin up density, now symmetrized
2047 49999 : call fourdp(cplex,rhog,work,1,mpi_enreg,nfft,1,ngfft,0)
2048 745284277 : rhor(:,ispden)=work(:)
2049 95253 : if (nspden==4) then
2050 1814 : call fourdp(cplex,magngx,work,1,mpi_enreg,nfft,1,ngfft,0)
2051 13723743 : rhor(:,2)=work(:)
2052 1814 : call fourdp(cplex,magngy,work,1,mpi_enreg,nfft,1,ngfft,0)
2053 13723743 : rhor(:,3)=work(:)
2054 1814 : call fourdp(cplex,magngz,work,1,mpi_enreg,nfft,1,ngfft,0)
2055 13723743 : rhor(:,4)=work(:)
2056 1814 : ABI_FREE(magngx)
2057 1814 : ABI_FREE(magngy)
2058 1814 : ABI_FREE(magngz)
2059 : end if
2060 :
2061 : end do ! ispden
2062 :
2063 : end if ! End on the condition nsym==1
2064 :
2065 142170 : ABI_FREE(work)
2066 :
2067 : contains
2068 :
2069 15357856 : integer function map_symrhg(j1, n1)
2070 :
2071 : integer,intent(in) :: j1, n1
2072 : ! Map into [0,n-1]
2073 15357856 : 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 6548 : 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 6548 : integer,allocatable :: class(:),iperm(:),symafm_used(:),symrel_used(:,:,:)
2153 6548 : integer,allocatable :: work1(:)
2154 6548 : real(dp),allocatable :: tnons_used(:,:),work2(:,:)
2155 : ! *************************************************************************
2156 :
2157 19644 : ABI_MALLOC(class,(nsym))
2158 13096 : ABI_MALLOC(iperm,(nsym))
2159 19644 : ABI_MALLOC(work1,(n1*n2*n3))
2160 19644 : ABI_MALLOC(work2,(2,n1*n2*n3))
2161 :
2162 6548 : nfftot=n1*n2*n3
2163 :
2164 6548 : id1=n1/2+2
2165 6548 : id2=n2/2+2
2166 6548 : id3=n3/2+2
2167 :
2168 6548 : sppoldbl=nspden/nsppol;if (nspden==4) sppoldbl=1
2169 :
2170 13127 : do imagn=1,sppoldbl
2171 :
2172 : ! Treat in a similar way the case of the full group and the non-magnetic subgroup
2173 6579 : nsym_used=0
2174 140951 : do isym=1,nsym
2175 134372 : if( (imagn==1 .and. sppoldbl==2) .or. symafm(isym)==1 .or. &
2176 6579 : & ((nspden==4).and.afm_noncoll) )then
2177 133556 : nsym_used=nsym_used+1
2178 : end if
2179 : end do
2180 :
2181 6579 : 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 19737 : ABI_MALLOC(symafm_used,(nsym_used))
2190 19737 : ABI_MALLOC(symrel_used,(3,3,nsym_used))
2191 19737 : ABI_MALLOC(tnons_used,(3,nsym_used))
2192 :
2193 6579 : nsym_used=0
2194 140951 : do isym=1,nsym
2195 134372 : if( (imagn==1 .and. sppoldbl==2) .or. symafm(isym)==1 .or. &
2196 6579 : & ((nspden==4).and.afm_noncoll) ) then
2197 133556 : nsym_used=nsym_used+1
2198 1736228 : symrel_used(:,:,nsym_used)=symrel(:,:,isym)
2199 534224 : tnons_used(:,nsym_used)=tnons(:,isym)
2200 133556 : symafm_used(nsym_used)=symafm(isym)
2201 : end if
2202 : end do
2203 139077 : 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 81829004 : work1(1:nfftot)=0
2209 81829004 : irrzon(:,2,imagn)=0
2210 :
2211 : ! Initialize at G=0 (always in irreducible zone)
2212 6579 : nzone=1
2213 6579 : irrzon(1,1,imagn)=1
2214 6579 : irrzon(1,2,imagn)=nsym_used
2215 : ! Set phase exp(2*Pi*I*G dot tnons) for G=0
2216 6579 : phnons(1,1,imagn)=one
2217 6579 : phnons(2,1,imagn)=zero
2218 6579 : npt=1
2219 : ! setting work1(1)=1 indicates that first grid point (G=0) is
2220 : ! in the iz (irreducible zone)
2221 6579 : work1(1)=1
2222 :
2223 6579 : ind1=0
2224 :
2225 : ! Loop over reciprocal space grid points:
2226 132199 : do i3=1,n3
2227 2856739 : do i2=1,n2
2228 84672585 : do i1=1,n1
2229 :
2230 81822425 : 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 84546965 : 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 10797138 : l3=i3-(i3/id3)*n3-1
2242 10797138 : l2=i2-(i2/id2)*n2-1
2243 10797138 : l1=i1-(i1/id1)*n1-1
2244 :
2245 116517022 : 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 105719884 : & symrel_used(2,1,isym)*l2+symrel_used(3,1,isym)*l3
2253 : j2=symrel_used(1,2,isym)*l1+&
2254 105719884 : & symrel_used(2,2,isym)*l2+symrel_used(3,2,isym)*l3
2255 : j3=symrel_used(1,3,isym)*l1+&
2256 105719884 : & symrel_used(2,3,isym)*l2+symrel_used(3,3,isym)*l3
2257 :
2258 : ! Map into [0,n-1]
2259 105719884 : k1=mod(n1+mod(j1,n1),n1)
2260 105719884 : k2=mod(n2+mod(j2,n2),n2)
2261 105719884 : k3=mod(n3+mod(j3,n3),n3)
2262 :
2263 : ! Get linear index of rotated point Gj
2264 105719884 : ind2=1+k1+n1*(k2+n2*k3)
2265 :
2266 : ! Store info for new class:
2267 105719884 : class(isym)=ind2
2268 105719884 : iperm(isym)=isym
2269 :
2270 : ! Setting work array element to 1 indicates grid point has been
2271 : ! identified with iz point
2272 116517022 : 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 10797138 : 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 10797138 : irep=0
2285 116517022 : do isym=1,nsym_used
2286 116517022 : if (class(isym)==class(1)) then
2287 14032170 : irep=irep+1
2288 : end if
2289 : end do
2290 10797138 : ipt=nsym_used/irep
2291 :
2292 : ! Repetition factor must be divisor of nsym_used:
2293 10797138 : 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 116517022 : 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 105719884 : isym=iperm(jj)
2310 105719884 : tau1=tnons_used(1,isym)
2311 105719884 : tau2=tnons_used(2,isym)
2312 105719884 : tau3=tnons_used(3,isym)
2313 : if (abs(tau1)>tol12.or.abs(tau2)>tol12&
2314 116517022 : & .or.abs(tau3)>tol12) then
2315 : ! compute exp(-2*Pi*I*G dot tau) using original G
2316 17053217 : arg=two_pi*(dble(l1)*tau1+dble(l2)*tau2+dble(l3)*tau3)
2317 17053217 : work2(1,jj)=cos(arg)
2318 17053217 : work2(2,jj)=-sin(arg)
2319 : else
2320 88666667 : work2(1,jj)=one
2321 88666667 : 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 10797138 : if (irep > 1) then
2346 17958482 : do jj=1,nsym_used,irep
2347 15748352 : setzer=0
2348 15748352 : ph1r=work2(1,jj);ph1i=work2(2,jj)
2349 55400742 : do j1=jj,jj+irep-1
2350 39652390 : ph2r=work2(1,j1);ph2i=work2(2,j1)
2351 55400742 : if (((ph2r+ph1r)**2+(ph2i+ph1i)**2) <= tol14) then
2352 39672582 : if (setzer/=1) setzer=-1
2353 36238508 : else if (((ph2r-ph1r)**2+(ph2i-ph1i)**2) > tol14) then
2354 196711 : 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 17958482 : if (setzer/=0) then
2361 1671992 : if (setzer==-1) then
2362 1571637 : if (afm_noncoll.and.nspden==4) then
2363 19484 : arg=symafm_used(iperm(jj))
2364 47879 : if (all(symafm_used(iperm(jj:jj+irep-1))==arg)) then
2365 : setzer=1
2366 : else
2367 40244 : do j1=jj,jj+irep-1
2368 94004 : 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 22910876 : if (setzer==1) work2(:,jj:jj+irep-1)=zero
2376 : end if
2377 : end do
2378 : ! Compress data if irep>1:
2379 2210130 : jj=0
2380 2210130 : do isym=1,nsym_used,irep
2381 15748352 : jj=jj+1
2382 15748352 : class(jj)=class(isym)
2383 15748352 : work2(1,jj)=work2(1,isym)
2384 15748352 : work2(2,jj)=work2(2,isym)
2385 : end do
2386 : end if
2387 :
2388 : ! Put new unique points into irrzon array:
2389 92612984 : irrzon(1+npt:ipt+npt,1,imagn)=class(1:ipt)
2390 :
2391 : ! Put repetition number into irrzon array:
2392 10797138 : 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 256244676 : 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 10797138 : npt=npt+ipt
2405 :
2406 : ! Update number of classes:
2407 10797138 : 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 6579 : ABI_SFREE(symafm_used)
2417 6579 : ABI_SFREE(symrel_used)
2418 13127 : 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 6548 : 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 13127 : do imagn=1,sppoldbl
2441 :
2442 81829004 : do ifft=1,nfftot
2443 81829004 : 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 6579 : izonemax=0
2452 81835552 : do izone=1,nfftot
2453 : ! Global bounds
2454 81822425 : 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 81822425 : if(izonemax==0)then
2462 10810296 : if (irrzon(izone,2,imagn)==0)izonemax=izone-1
2463 : end if
2464 10816875 : if(izonemax/=0)then
2465 71018708 : 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 6548 : ABI_FREE(class)
2475 6548 : ABI_FREE(iperm)
2476 6548 : ABI_FREE(work1)
2477 6548 : ABI_FREE(work2)
2478 :
2479 6548 : 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 9423 : subroutine setsym(indsym,irrzon,iscf,natom,nfft,ngfft,nspden,nsppol,nsym,phnons,&
2681 9423 : & 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 9423 : integer,allocatable :: determinant(:)
2700 : real(dp) :: tsec(2)
2701 : ! *************************************************************************
2702 :
2703 9423 : call timab(6,1,tsec)
2704 :
2705 : !Check that symmetries have unity determinant
2706 28269 : ABI_MALLOC(determinant,(nsym))
2707 9423 : call symdet(determinant,nsym,symrel)
2708 9423 : ABI_FREE(determinant)
2709 :
2710 :
2711 : !Get the symmetry matrices in terms of reciprocal basis
2712 144701 : do isym=1,nsym
2713 144701 : call mati3inv(symrel(:,:,isym),symrec(:,:,isym))
2714 : end do
2715 :
2716 : !Check for group closure
2717 9423 : call sg_multable(nsym,symafm,symrel,ierr,tnons=tnons,tnons_tol=tol5)
2718 9423 : ABI_CHECK(ierr==0,"Error in group closure")
2719 :
2720 9423 : call sg_multable(nsym,symafm,symrec,ierr)
2721 9423 : ABI_CHECK(ierr==0,"Error in group closure")
2722 :
2723 : !Obtain a list of rotated atom labels:
2724 9423 : tolsym8=tol8
2725 9423 : 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 9423 : if ( (iscf>0 .or. iscf==-3) .and. nsym>1 ) then
2729 : ! Locate irreducible zone in reciprocal space for symmetrization:
2730 6089 : call irrzg(irrzon,nspden,nsppol,nsym,ngfft(1),ngfft(2),ngfft(3),phnons,symafm,symrel,tnons)
2731 : end if
2732 :
2733 9423 : call timab(6,2,tsec)
2734 :
2735 : !DEBUG
2736 : !write(std_out,*)' m_spacepar%setsym : exit '
2737 : !ENDDEBUG
2738 :
2739 9423 : 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 : !!***
|
