TABLE OF CONTENTS


ABINIT/m_spacepar [ Modules ]

[ Top ] [ Modules ]

NAME

 m_spacepar

FUNCTION

  Relatively Low-level procedures operating on arrays defined on the FFT box (G- or R- space)
  Unlike the procedures in m_cgtools, the routines declared in this module can use mpi_type.

COPYRIGHT

  Copyright (C) 2008-2018 ABINIT group (XG, BA, MT, DRH, DCA, GMR, MJV)
  This file is distributed under the terms of the
  GNU General Public License, see ~abinit/COPYING
  or http://www.gnu.org/copyleft/gpl.txt .

PARENTS

CHILDREN

SOURCE

22 #if defined HAVE_CONFIG_H
23 #include "config.h"
24 #endif
25 
26 #include "abi_common.h"
27 
28 module m_spacepar
29 
30  use defs_basis
31  use m_abicore
32  use m_errors
33  use m_xmpi
34  use m_sort
35 
36  use m_time,            only : timab
37  use defs_abitypes,     only : MPI_type
38  use m_symtk,           only : mati3inv, chkgrp, symdet, symatm, matr3inv
39  use m_geometry,        only : metric, symredcart
40  use m_mpinfo,          only : ptabs_fourdp
41  use m_fft,             only : zerosym, fourdp
42 
43  implicit none
44 
45  private

m_spacepar/hartre [ Functions ]

[ Top ] [ m_spacepar ] [ Functions ]

NAME

 hartre

FUNCTION

 Given rho(G), compute Hartree potential (=FFT of rho(G)/pi/(G+q)**2)
 When cplex=1, assume q=(0 0 0), and vhartr will be REAL
 When cplex=2, q must be taken into account, and vhartr will be COMPLEX

NOTES

 *Modified code to avoid if statements inside loops to skip G=0.
  Replaced if statement on G^2>gsqcut to skip G s outside where
  rho(G) should be 0.  Effect is negligible but gsqcut should be
  used to be strictly consistent with usage elsewhere in code.
 *The speed-up is provided by doing a few precomputations outside
  the inner loop. One variable size array is needed for this (gq).

INPUTS

  cplex= if 1, vhartr is REAL, if 2, vhartr is COMPLEX
  gsqcut=cutoff value on G**2 for sphere inside fft box.
         (gsqcut=(boxcut**2)*ecut/(2.d0*(Pi**2))
  izero=if 1, unbalanced components of Vhartree(g) are set to zero
  mpi_enreg=information about MPI parallelization
  nfft=(effective) number of FFT grid points (for this processor)
  ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft
  [qpt(3)=reduced coordinates for a wavevector to be combined with the G vectors (needed if cplex==2).]
  rhog(2,nfft)=electron density in G space
  rprimd(3,3)=dimensional primitive translations in real space (bohr)
  divgq0= [optional argument] value of the integration of the Coulomb singularity 4pi\int_BZ 1/q^2 dq

OUTPUT

  vhartr(cplex*nfft)=Hartree potential in real space, either REAL or COMPLEX

PARENTS

      dfpt_rhotov,dfptnl_loop,energy,fock_getghc,m_kxc,nonlinear,nres2vres
      odamix,prcref,prcref_PMA,respfn,rhotov,setup_positron,setvtr,tddft

CHILDREN

      fourdp,metric,ptabs_fourdp,timab,zerosym

SOURCE

106 subroutine hartre(cplex,gsqcut,izero,mpi_enreg,nfft,ngfft,paral_kgb,rhog,rprimd,vhartr,&
107 &  divgq0,qpt) ! Optional argument
108 
109 
110 !This section has been created automatically by the script Abilint (TD).
111 !Do not modify the following lines by hand.
112 #undef ABI_FUNC
113 #define ABI_FUNC 'hartre'
114 !End of the abilint section
115 
116  implicit none
117 
118 !Arguments ------------------------------------
119 !scalars
120  integer,intent(in) :: cplex,izero,nfft,paral_kgb
121  real(dp),intent(in) :: gsqcut
122  type(MPI_type),intent(in) :: mpi_enreg
123 !arrays
124  integer,intent(in) :: ngfft(18)
125  real(dp),intent(in) :: rprimd(3,3),rhog(2,nfft)
126  real(dp),intent(in),optional :: divgq0
127  real(dp),intent(in),optional :: qpt(3)
128  real(dp),intent(out) :: vhartr(cplex*nfft)
129 
130 !Local variables-------------------------------
131 !scalars
132  integer,parameter :: im=2,re=1
133  integer :: i1,i2,i23,i2_local,i3,id1,id2,id3
134  integer :: ig,ig1min,ig1,ig1max,ig2,ig2min,ig2max,ig3,ig3min,ig3max
135  integer :: ii,ii1,ing,n1,n2,n3,qeq0,qeq05,me_fft,nproc_fft
136  real(dp),parameter :: tolfix=1.000000001e0_dp
137  real(dp) :: cutoff,den,gqg2p3,gqgm12,gqgm13,gqgm23,gs,gs2,gs3,ucvol
138  character(len=500) :: message
139 !arrays
140  integer :: id(3)
141  integer, ABI_CONTIGUOUS pointer :: fftn2_distrib(:),ffti2_local(:)
142  integer, ABI_CONTIGUOUS pointer :: fftn3_distrib(:),ffti3_local(:)
143  real(dp) :: gmet(3,3),gprimd(3,3),qpt_(3),rmet(3,3),tsec(2)
144  real(dp),allocatable :: gq(:,:),work1(:,:)
145 
146 ! *************************************************************************
147 
148  ! Keep track of total time spent in hartre
149  call timab(10,1,tsec)
150 
151  ! Check that cplex has an allowed value
152  if(cplex/=1 .and. cplex/=2)then
153    write(message, '(a,i0,a,a)' )&
154 &   'From the calling routine, cplex=',cplex,ch10,&
155 &   'but the only value allowed are 1 and 2.'
156    MSG_BUG(message)
157  end if
158 
159  call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
160 
161  n1=ngfft(1); n2=ngfft(2); n3=ngfft(3)
162  nproc_fft = mpi_enreg%nproc_fft; me_fft = mpi_enreg%me_fft
163 
164  ! Get the distrib associated with this fft_grid
165  call ptabs_fourdp(mpi_enreg,n2,n3,fftn2_distrib,ffti2_local,fftn3_distrib,ffti3_local)
166 
167  ! Initialize a few quantities
168  cutoff=gsqcut*tolfix
169  if(present(qpt))then
170    qpt_=qpt
171  else
172    qpt_=zero
173  end if
174  qeq0=0
175  if(qpt_(1)**2+qpt_(2)**2+qpt_(3)**2<1.d-15) qeq0=1
176  qeq05=0
177  if (qeq0==0) then
178    if (abs(abs(qpt_(1))-half)<tol12.or.abs(abs(qpt_(2))-half)<tol12.or. &
179 &   abs(abs(qpt_(3))-half)<tol12) qeq05=1
180  end if
181 
182  ! If cplex=1 then qpt_ should be 0 0 0
183  if (cplex==1.and. qeq0/=1) then
184    write(message,'(a,3e12.4,a,a)')&
185 &   'cplex=1 but qpt=',qpt_,ch10,&
186 &   'qpt should be 0 0 0.'
187    MSG_BUG(message)
188  end if
189 
190  ! If FFT parallelism then qpt should not be 1/2
191  if (nproc_fft>1.and.qeq05==1) then
192    write(message, '(a,3e12.4,a,a)' )&
193 &   'FFT parallelism selected but qpt',qpt_,ch10,&
194 &   'qpt(i) should not be 1/2...'
195    MSG_ERROR(message)
196  end if
197 
198  ! In order to speed the routine, precompute the components of g+q
199  ! Also check if the booked space was large enough...
200  ABI_ALLOCATE(gq,(3,max(n1,n2,n3)))
201  do ii=1,3
202    id(ii)=ngfft(ii)/2+2
203    do ing=1,ngfft(ii)
204      ig=ing-(ing/id(ii))*ngfft(ii)-1
205      gq(ii,ing)=ig+qpt_(ii)
206    end do
207  end do
208  ig1max=-1;ig2max=-1;ig3max=-1
209  ig1min=n1;ig2min=n2;ig3min=n3
210 
211  ABI_ALLOCATE(work1,(2,nfft))
212  id1=n1/2+2;id2=n2/2+2;id3=n3/2+2
213 
214  ! Triple loop on each dimension
215  do i3=1,n3
216    ig3=i3-(i3/id3)*n3-1
217    ! Precompute some products that do not depend on i2 and i1
218    gs3=gq(3,i3)*gq(3,i3)*gmet(3,3)
219    gqgm23=gq(3,i3)*gmet(2,3)*2
220    gqgm13=gq(3,i3)*gmet(1,3)*2
221 
222    do i2=1,n2
223      ig2=i2-(i2/id2)*n2-1
224      if (fftn2_distrib(i2) == me_fft) then
225        gs2=gs3+ gq(2,i2)*(gq(2,i2)*gmet(2,2)+gqgm23)
226        gqgm12=gq(2,i2)*gmet(1,2)*2
227        gqg2p3=gqgm13+gqgm12
228 
229        i2_local = ffti2_local(i2)
230        i23=n1*(i2_local-1 +(n2/nproc_fft)*(i3-1))
231        ! Do the test that eliminates the Gamma point outside of the inner loop
232        ii1=1
233        if(i23==0 .and. qeq0==1  .and. ig2==0 .and. ig3==0)then
234          ii1=2
235          work1(re,1+i23)=zero
236          work1(im,1+i23)=zero
237          ! If the value of the integration of the Coulomb singularity 4pi\int_BZ 1/q^2 dq is given, use it
238          if (PRESENT(divgq0)) then
239            work1(re,1+i23)=rhog(re,1+i23)*divgq0*piinv
240            work1(im,1+i23)=rhog(im,1+i23)*divgq0*piinv
241          end if
242        end if
243 
244        ! Final inner loop on the first dimension (note the lower limit)
245        do i1=ii1,n1
246          gs=gs2+ gq(1,i1)*(gq(1,i1)*gmet(1,1)+gqg2p3)
247          ii=i1+i23
248 
249          if(gs<=cutoff)then
250            ! Identify min/max indexes (to cancel unbalanced contributions later)
251            ! Count (q+g)-vectors with similar norm
252            if ((qeq05==1).and.(izero==1)) then
253              ig1=i1-(i1/id1)*n1-1
254              ig1max=max(ig1max,ig1); ig1min=min(ig1min,ig1)
255              ig2max=max(ig2max,ig2); ig2min=min(ig2min,ig2)
256              ig3max=max(ig3max,ig3); ig3min=min(ig3min,ig3)
257            end if
258 
259            den=piinv/gs
260            work1(re,ii)=rhog(re,ii)*den
261            work1(im,ii)=rhog(im,ii)*den
262          else
263            ! gs>cutoff
264            work1(re,ii)=zero
265            work1(im,ii)=zero
266          end if
267 
268        end do ! End loop on i1
269      end if
270    end do ! End loop on i2
271  end do ! End loop on i3
272 
273  ABI_DEALLOCATE(gq)
274 
275  if (izero==1) then
276    ! Set contribution of unbalanced components to zero
277 
278    if (qeq0==1) then !q=0
279      call zerosym(work1,2,n1,n2,n3,comm_fft=mpi_enreg%comm_fft,distribfft=mpi_enreg%distribfft)
280 
281    else if (qeq05==1) then
282      !q=1/2; this doesn't work in parallel
283      ig1=-1;if (mod(n1,2)==0) ig1=1+n1/2
284      ig2=-1;if (mod(n2,2)==0) ig2=1+n2/2
285      ig3=-1;if (mod(n3,2)==0) ig3=1+n3/2
286      if (abs(abs(qpt_(1))-half)<tol12) then
287        if (abs(ig1min)<abs(ig1max)) ig1=abs(ig1max)
288        if (abs(ig1min)>abs(ig1max)) ig1=n1-abs(ig1min)
289      end if
290      if (abs(abs(qpt_(2))-half)<tol12) then
291        if (abs(ig2min)<abs(ig2max)) ig2=abs(ig2max)
292        if (abs(ig2min)>abs(ig2max)) ig2=n2-abs(ig2min)
293      end if
294      if (abs(abs(qpt_(3))-half)<tol12) then
295        if (abs(ig3min)<abs(ig3max)) ig3=abs(ig3max)
296        if (abs(ig3min)>abs(ig3max)) ig3=n3-abs(ig3min)
297      end if
298      call zerosym(work1,2,n1,n2,n3,ig1=ig1,ig2=ig2,ig3=ig3,&
299 &     comm_fft=mpi_enreg%comm_fft,distribfft=mpi_enreg%distribfft)
300    end if
301  end if
302 
303  ! Fourier Transform Vhartree. Vh in reciprocal space was stored in work1
304  call fourdp(cplex,work1,vhartr,1,mpi_enreg,nfft,ngfft,paral_kgb,0)
305 
306  ABI_DEALLOCATE(work1)
307 
308  call timab(10,2,tsec)
309 
310 end subroutine hartre

m_spacepar/hartrestr [ Functions ]

[ Top ] [ m_spacepar ] [ Functions ]

NAME

 hartrestr

FUNCTION

 To be called for strain perturbation only
 Compute the inhomogenous terms generated by the strain derivative of
 Hartree potential due to the ground state charge rho(G)

  FFT of (rho(G)/pi)*[d(1/G**2)/d(strain) - delta(diagonal strain)*(1/G**2)]

NOTES

 *based largely on hartre.f
 *Modified code to avoid if statements inside loops to skip G=0.
  Replaced if statement on G^2>gsqcut to skip G s outside where
  rho(G) should be 0.  Effect is negligible but gsqcut should be
  used to be strictly consistent with usage elsewhere in code.
 *The speed-up is provided by doing a few precomputations outside
  the inner loop. One variable size array is needed for this (gq).

INPUTS

  gsqcut=cutoff value on G**2 for sphere inside fft box.
  idir=direction of the current perturbation
  ipert=type of the perturbation
  mpi_enreg=information about MPI parallelization
  natom=number of atoms in cell.
  nfft=number of fft grid points (gsqcut=(boxcut**2)*ecut/(2._dp*(Pi**2))
  ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft
  rhog(2,nfft)=array for Fourier transform of GS electron density
  rprimd(3,3)=dimensional primitive translations in real space (bohr)

OUTPUT

  vhartr1(nfft)=Inhomogeneous term in strain-perturbation-induced Hartree
   potential in real space,

PARENTS

      dfpt_nselt,dfpt_nstpaw,dfpt_rhotov

CHILDREN

      fourdp,metric,ptabs_fourdp

SOURCE

 911 subroutine hartrestr(gsqcut,idir,ipert,mpi_enreg,natom,nfft,ngfft,&
 912 &  paral_kgb,rhog,rprimd,vhartr1)
 913 
 914 
 915 !This section has been created automatically by the script Abilint (TD).
 916 !Do not modify the following lines by hand.
 917 #undef ABI_FUNC
 918 #define ABI_FUNC 'hartrestr'
 919 !End of the abilint section
 920 
 921  implicit none
 922 
 923 !Arguments ------------------------------------
 924 !scalars
 925  integer,intent(in) :: idir,ipert,natom,nfft,paral_kgb
 926  real(dp),intent(in) :: gsqcut
 927  type(MPI_type),intent(in) :: mpi_enreg
 928 !arrays
 929  integer,intent(in) :: ngfft(18)
 930  real(dp),intent(in) :: rhog(2,nfft),rprimd(3,3)
 931  real(dp),intent(out) :: vhartr1(nfft)
 932 
 933 !Local variables-------------------------------
 934 !scalars
 935  integer,parameter :: im=2,re=1
 936  integer :: i1,i2,i23,i3,id2,id3,ig,ig2,ig3,ii,ii1,ing,istr,ka,kb,n1,n2,n3
 937  real(dp),parameter :: tolfix=1.000000001_dp
 938  real(dp) :: cutoff,ddends,den,dgsds,gqg2p3,gqgm12,gqgm13,gqgm23,gs,gs2,gs3
 939  real(dp) :: term,ucvol
 940  character(len=500) :: message
 941 !arrays
 942  integer,save :: idx(12)=(/1,1,2,2,3,3,3,2,3,1,2,1/)
 943  integer :: id(3)
 944  integer, ABI_CONTIGUOUS pointer :: fftn2_distrib(:),ffti2_local(:)
 945  integer, ABI_CONTIGUOUS pointer :: fftn3_distrib(:),ffti3_local(:)
 946  real(dp) :: dgmetds(3,3),gmet(3,3),gprimd(3,3),gqr(3),rmet(3,3)
 947  real(dp),allocatable :: gq(:,:),work1(:,:)
 948 
 949 ! *************************************************************************
 950 
 951  if( .not. (ipert==natom+3 .or. ipert==natom+4))then
 952    write(message, '(a,i0,a,a)' )&
 953 &   'From the calling routine, ipert=',ipert,ch10,&
 954 &   'so this routine for the strain perturbation should not be called.'
 955    MSG_BUG(message)
 956  end if
 957 
 958  call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
 959 
 960  n1=ngfft(1) ; n2=ngfft(2) ; n3=ngfft(3)
 961 
 962 !Get the distrib associated with this fft_grid
 963  call ptabs_fourdp(mpi_enreg,n2,n3,fftn2_distrib,ffti2_local,fftn3_distrib,ffti3_local)
 964 
 965 !Initialize a few quantities
 966  cutoff=gsqcut*tolfix
 967 
 968  istr=idir + 3*(ipert-natom-3)
 969 
 970  if(istr<1 .or. istr>6)then
 971    write(message, '(a,i10,a,a,a)' )&
 972 &   'Input dir gives istr=',istr,' not allowed.',ch10,&
 973 &   'Possible values are 1,2,3,4,5,6 only.'
 974    MSG_BUG(message)
 975  end if
 976 
 977  ka=idx(2*istr-1);kb=idx(2*istr)
 978  do ii = 1,3
 979    dgmetds(:,ii)=-(gprimd(ka,:)*gprimd(kb,ii)+gprimd(kb,:)*gprimd(ka,ii))
 980  end do
 981 !For historical reasons:
 982  dgmetds(:,:)=0.5_dp*dgmetds(:,:)
 983 
 984 !In order to speed the routine, precompute the components of g+q
 985 !Also check if the booked space was large enough...
 986  ABI_ALLOCATE(gq,(3,max(n1,n2,n3)))
 987  do ii=1,3
 988    id(ii)=ngfft(ii)/2+2
 989    do ing=1,ngfft(ii)
 990      ig=ing-(ing/id(ii))*ngfft(ii)-1
 991      gq(ii,ing)=ig
 992    end do
 993  end do
 994 
 995  ABI_ALLOCATE(work1,(2,nfft))
 996  id2=n2/2+2
 997  id3=n3/2+2
 998 !Triple loop on each dimension
 999  do i3=1,n3
1000    ig3=i3-(i3/id3)*n3-1
1001 !  Precompute some products that do not depend on i2 and i1
1002    gqr(3)=gq(3,i3)
1003    gs3=gq(3,i3)*gq(3,i3)*gmet(3,3)
1004    gqgm23=gq(3,i3)*gmet(2,3)*2
1005    gqgm13=gq(3,i3)*gmet(1,3)*2
1006 
1007    do i2=1,n2
1008      if (fftn2_distrib(i2)==mpi_enreg%me_fft) then
1009        gqr(2)=gq(2,i2)
1010        gs2=gs3+ gq(2,i2)*(gq(2,i2)*gmet(2,2)+gqgm23)
1011        gqgm12=gq(2,i2)*gmet(1,2)*2
1012        gqg2p3=gqgm13+gqgm12
1013        ig2=i2-(i2/id2)*n2-1
1014 !      i23=n1*((i2-1)+n2*(i3-1))
1015        i23=n1*((ffti2_local(i2)-1)+(n2/mpi_enreg%nproc_fft)*(i3-1))
1016 !      Do the test that eliminates the Gamma point outside
1017 !      of the inner loop
1018        ii1=1
1019        if(i23==0  .and. ig2==0 .and. ig3==0)then
1020          ii1=2
1021          work1(re,1+i23)=0.0_dp
1022          work1(im,1+i23)=0.0_dp
1023        end if
1024 
1025 !      Final inner loop on the first dimension
1026 !      (note the lower limit)
1027        do i1=ii1,n1
1028          gs=gs2+ gq(1,i1)*(gq(1,i1)*gmet(1,1)+gqg2p3)
1029          ii=i1+i23
1030          if(gs<=cutoff)then
1031            den=piinv/gs
1032            gqr(1)=gq(1,i1)
1033            dgsds=&
1034 &           (gqr(1)*(dgmetds(1,1)*gqr(1)+dgmetds(1,2)*gqr(2)+dgmetds(1,3)*gqr(3))+  &
1035 &           gqr(2)*(dgmetds(2,1)*gqr(1)+dgmetds(2,2)*gqr(2)+dgmetds(2,3)*gqr(3))+  &
1036 &           gqr(3)*(dgmetds(3,1)*gqr(1)+dgmetds(3,2)*gqr(2)+dgmetds(3,3)*gqr(3)) )
1037            ddends=-piinv*dgsds/gs**2
1038            if(istr<=3)then
1039              term=2.0_dp*ddends-den
1040            else
1041              term=2.0_dp*ddends
1042            end if
1043            work1(re,ii)=rhog(re,ii)*term
1044            work1(im,ii)=rhog(im,ii)*term
1045          else
1046            work1(re,ii)=0.0_dp
1047            work1(im,ii)=0.0_dp
1048          end if
1049 
1050        end do ! End loop on i1
1051      end if
1052    end do ! End loop on i2
1053  end do !  End loop on i3
1054 
1055  ABI_DEALLOCATE(gq)
1056 
1057 !Fourier Transform Vhartree.
1058 !Vh in reciprocal space was stored in work1
1059  call fourdp(1,work1,vhartr1,1,mpi_enreg,nfft,ngfft,paral_kgb,0)
1060 
1061  ABI_DEALLOCATE(work1)
1062 
1063 end subroutine hartrestr

m_spacepar/irrzg [ Functions ]

[ Top ] [ m_spacepar ] [ Functions ]

NAME

 irrzg

FUNCTION

 Find the irreducible zone in reciprocal space under the
 symmetry group with real space rotations in symrel(3,3,nsym).
 The (integer) rotation matrices symrel(3,3,nsym) express the new
 real space positions (e.g. rotated atom positions) in REDUCED
 coordinates, i.e. in coordinates expressed as fractions of real space
 primitive translations (atomic coordinates xred).  tnons(3,nsym) express
 the associated nonsymmorphic translations, again in reduced coordinates.
 Special data structure created in irrzon.
 First half holds mapping from irr zone to full zone;
 part of second half holds repetition number info.
 work1 is a work array to keep track of grid points found so far.
 In case nspden=2 and nsppol=1, one has to take care of antiferromagnetic
 operations. The subgroup of non-magnetic operations is used
 to generate irrzon(:,:,2) and phnons(:,:,2), while the
 full group is used to generate irrzon(:,:,1) and phnons(:,:,1)

NOTES

 for reference in the near future (2018), some notes: this routine should be duplicated for
 magnetizations in spinorial formalism. The only main difference will
 be that the density is not simply transported to the image point under symrel
 but is a vector which has to be transformed as well
 $S \vec{m} (\vec{x}) = \sigma \vec{m} (S \vec{x} + \tau)$
 $S \vec{m} (\vec{G}) = \sigma exp(+ 2 \pi i \vec{G} \vec{tau}) \vec{m}(S^{-1 T} \vec{G})$
 S is a symop, sigma the AFM sign flip if any, tau the partial non symmorphic translation
 x a position, m the magnetization 3 vector

 The phase factor is the same as below for the density, but the collection of elements which
 are equal is more complex: for a 3 or 6 axis the m vector could transform one component
 to a linear combination of several others (I think). Things are not necessarily aligned
 with the z axis.

INPUTS

  nspden=number of spin-density components
  nsppol=1 for unpolarized, 2 for spin-polarized
  nsym=number of symmetry elements in group
  n1,n2,n3=box dimensions of real space grid (or fft box)
  symafm(nsym)=(anti)ferromagnetic part of symmetry operations
  symrel(3,3,nsym)=symmetry matrices in real space (integers)
  tnons(3,nsym)=reduced nonsymmorphic translations
 (symrel and tnons are in terms of real space primitive translations)

OUTPUT

  irrzon(n1*n2*n3,2+(nspden/4),(nspden/nsppol)-3*(nspden/4))=integer array which contains the locations of related
   grid points and the repetition number for each symmetry class.
  phnons(2,n1*n2*n3,(nspden/nsppol)-3*(nspden/4))=phases associated with nonsymmorphic translations

PARENTS

      m_ab7_kpoints,setsym,wfd_mkrho

CHILDREN

      sort_int,wrtout

SOURCE

1575 subroutine irrzg(irrzon,nspden,nsppol,nsym,n1,n2,n3,phnons,symafm,symrel,tnons)
1576 
1577 
1578 !This section has been created automatically by the script Abilint (TD).
1579 !Do not modify the following lines by hand.
1580 #undef ABI_FUNC
1581 #define ABI_FUNC 'irrzg'
1582 !End of the abilint section
1583 
1584  implicit none
1585 
1586 !Arguments ------------------------------------
1587 !scalars
1588  integer,intent(in) :: n1,n2,n3,nspden,nsppol,nsym
1589 !arrays
1590  integer,intent(in) :: symafm(nsym),symrel(3,3,nsym)
1591  integer,intent(out) :: irrzon(n1*n2*n3,2,(nspden/nsppol)-3*(nspden/4))
1592  real(dp),intent(in) :: tnons(3,nsym)
1593  real(dp),intent(out) :: phnons(2,n1*n2*n3,(nspden/nsppol)-3*(nspden/4))
1594 
1595 !Local variables-------------------------------
1596 !scalars
1597  integer :: i1,i2,i3,id1,id2,id3,ifft,imagn,ind1,ind2,ipt,irep,isym,izone
1598  integer :: izonemax,j1,j2,j3,jj,k1,k2,k3,l1,l2,l3,nfftot,npt,nsym_used
1599  integer :: nzone,setzer,sppoldbl
1600  real(dp) :: arg,ph1i,ph1r,ph2i,ph2r,tau1,tau2,tau3
1601  logical,parameter :: afm_noncoll=.true. ! TRUE if antiferro symmetries are used in non-collinear magnetism
1602  character(len=500) :: message
1603 !arrays
1604  integer,allocatable :: class(:),iperm(:),symafm_used(:),symrel_used(:,:,:)
1605  integer,allocatable :: work1(:)
1606  real(dp),allocatable :: tnons_used(:,:),work2(:,:)
1607 
1608 ! *************************************************************************
1609 
1610  ABI_ALLOCATE(class,(nsym))
1611  ABI_ALLOCATE(iperm,(nsym))
1612  ABI_ALLOCATE(work1,(n1*n2*n3))
1613  ABI_ALLOCATE(work2,(2,n1*n2*n3))
1614 
1615  nfftot=n1*n2*n3
1616 
1617  id1=n1/2+2
1618  id2=n2/2+2
1619  id3=n3/2+2
1620 
1621  sppoldbl=nspden/nsppol;if (nspden==4) sppoldbl=1
1622 
1623  do imagn=1,sppoldbl
1624 
1625 !  Treat in a similar way the case of the full group and the non-magnetic subgroup
1626    nsym_used=0
1627    do isym=1,nsym
1628      if( (imagn==1 .and. sppoldbl==2) .or. symafm(isym)==1 .or. &
1629 &     ((nspden==4).and.afm_noncoll) )then
1630        nsym_used=nsym_used+1
1631      end if
1632    end do
1633 
1634    if(imagn==2 .and. nsym_used/=nsym/2)then
1635      write(message, '(a,a,a,a,a,i4,a,i0)' )&
1636 &     '  The number of ferromagnetic symmetry operations must be',ch10,&
1637 &     '  half the total number of operations, while it is observed that',ch10,&
1638 &     '  nsym=',nsym,' and nsym_magn=',nsym_used
1639      MSG_BUG(message)
1640    end if
1641 
1642    ABI_ALLOCATE(symafm_used,(nsym_used))
1643    ABI_ALLOCATE(symrel_used,(3,3,nsym_used))
1644    ABI_ALLOCATE(tnons_used,(3,nsym_used))
1645 
1646    nsym_used=0
1647    do isym=1,nsym
1648      if( (imagn==1 .and. sppoldbl==2) .or. symafm(isym)==1 .or.  &
1649 &     ((nspden==4).and.afm_noncoll) ) then
1650        nsym_used=nsym_used+1
1651        symrel_used(:,:,nsym_used)=symrel(:,:,isym)
1652        tnons_used(:,nsym_used)=tnons(:,isym)
1653        symafm_used(nsym_used)=symafm(isym)
1654      end if
1655    end do
1656    if ((nspden/=4).or.(.not.afm_noncoll)) symafm_used=1
1657 
1658 
1659 !  Zero out work array--later on, a zero entry will mean that
1660 !  a given grid point has not yet been assigned to an ibz point
1661    work1(1:nfftot)=0
1662    irrzon(:,2,imagn)=0
1663 
1664 !  Initialize at G=0 (always in irreducible zone)
1665    nzone=1
1666    irrzon(1,1,imagn)=1
1667    irrzon(1,2,imagn)=nsym_used
1668 !  Set phase exp(2*Pi*I*G dot tnons) for G=0
1669    phnons(1,1,imagn)=one
1670    phnons(2,1,imagn)=zero
1671    npt=1
1672 !  setting work1(1)=1 indicates that first grid point (G=0) is
1673 !  in the iz (irreducible zone)
1674    work1(1)=1
1675 
1676    ind1=0
1677 
1678 !  Loop over reciprocal space grid points:
1679    do i3=1,n3
1680      do i2=1,n2
1681        do i1=1,n1
1682 
1683          ind1=ind1+1
1684 
1685 !        Check to see whether present grid point is equivalent to
1686 !        any previously identified ibz point--if not, a new ibz point
1687 !        has been found
1688 
1689          if (work1(ind1)==0) then
1690 
1691 !          A new class has been found.
1692 
1693 !          Get location of G vector (grid point) centered at 0 0 0
1694            l3=i3-(i3/id3)*n3-1
1695            l2=i2-(i2/id2)*n2-1
1696            l1=i1-(i1/id1)*n1-1
1697 
1698            do isym=1,nsym_used
1699 
1700 !            Get rotated G vector Gj for each symmetry element
1701 !            -- here we use the TRANSPOSE of symrel_used; assuming symrel_used expresses
1702 !            the rotation in real space, the transpose is then appropriate
1703 !            for G space symmetrization (p. 1172d,e of notes, 2 June 1995).
1704              j1=symrel_used(1,1,isym)*l1+&
1705 &             symrel_used(2,1,isym)*l2+symrel_used(3,1,isym)*l3
1706              j2=symrel_used(1,2,isym)*l1+&
1707 &             symrel_used(2,2,isym)*l2+symrel_used(3,2,isym)*l3
1708              j3=symrel_used(1,3,isym)*l1+&
1709 &             symrel_used(2,3,isym)*l2+symrel_used(3,3,isym)*l3
1710 
1711 !            Map into [0,n-1] and then add 1 for array index in [1,n]
1712              k1=1+mod(n1+mod(j1,n1),n1)
1713              k2=1+mod(n2+mod(j2,n2),n2)
1714              k3=1+mod(n3+mod(j3,n3),n3)
1715 !            k1=1+map(j1,n1)
1716 !            k2=1+map(j2,n2)
1717 !            k3=1+map(j3,n3)
1718 
1719 !            Get linear index of rotated point Gj
1720              ind2=k1+n1*((k2-1)+n2*(k3-1))
1721 
1722 !            Store info for new class:
1723              class(isym)=ind2
1724              iperm(isym)=isym
1725 
1726 !            Setting work array element to 1 indicates grid point has been
1727 !            identified with iz point
1728              work1(ind2)=1
1729 
1730 !            End of loop on isym
1731            end do
1732 
1733 !          Sort integers into ascending order in each class
1734 !          (this lumps together G vectors with the same linear index, i.e.
1735 !          groups together symmetries which land on the same Gj)
1736            call sort_int(nsym_used,class,iperm)
1737 
1738 !          Check repetition factor (how many identical copies of Gj occur
1739 !          from all symmetries applied to G)
1740            irep=0
1741            do isym=1,nsym_used
1742              if (class(isym)==class(1)) then
1743                irep=irep+1
1744              end if
1745            end do
1746            ipt=nsym_used/irep
1747 
1748 !          Repetition factor must be divisor of nsym_used:
1749            if (nsym_used/=(ipt*irep)) then
1750              write(message, '(a,i5,a,i6,a,a,a,a,a,a)' )&
1751 &             '  irep=',irep,' not a divisor of nsym_used=',nsym_used,ch10,&
1752 &             ' This usually indicates that',&
1753 &             ' the input symmetries do not form a group.',ch10,&
1754 &             ' Action : check the input symmetries carefully do they',&
1755 &             ' form a group ? If they do, there is a code bug.'
1756              MSG_ERROR(message)
1757            end if
1758 
1759 !          Compute phases for any nonsymmorphic symmetries
1760 !          exp(-2*Pi*I*G dot tau(j)) for each symmetry j with
1761 !          (possibly zero) nonsymmorphic translation tau(j)
1762            do jj=1,nsym_used
1763 !            First get nonsymmorphic translation and see if nonzero
1764 !            (iperm grabs the symmetries in the new order after sorting)
1765              isym=iperm(jj)
1766              tau1=tnons_used(1,isym)
1767              tau2=tnons_used(2,isym)
1768              tau3=tnons_used(3,isym)
1769              if (abs(tau1)>tol12.or.abs(tau2)>tol12&
1770 &             .or.abs(tau3)>tol12) then
1771 !              compute exp(-2*Pi*I*G dot tau) using original G
1772                arg=two_pi*(dble(l1)*tau1+dble(l2)*tau2+dble(l3)*tau3)
1773                work2(1,jj)=cos(arg)
1774                work2(2,jj)=-sin(arg)
1775              else
1776                work2(1,jj)=one
1777                work2(2,jj)=zero
1778              end if
1779            end do
1780 
1781 !          All phases arising from symmetries which map to the same
1782 !          G vector must actually be the same because
1783 !          rho(Strans*G)=exp(2*Pi*I*(G) dot tau_S) rho(G)
1784 !          must be satisfied; if exp(2*Pi*I*(G) dot tau_S) can be different
1785 !          for two different symmetries S which both take G to the same St*G,
1786 !          then the related Fourier components rho(St*G) must VANISH.
1787 !          Thus: set "phase" to ZERO here in that case.
1788 !          The G mappings occur in sets of irep members; if irep=1 then
1789 !          all the G are unique.
1790 !          MT 090212:
1791 !          In the case of antiferromagn. symetries, one can have
1792 !          rho(Strans*G)= -exp(2*Pi*I*(G) dot tau_S) rho(G)
1793 !          (look at the minus !)
1794 !          A special treatment is then operated on phons.
1795 !          The later must be consistent with the use of phnons array
1796 !          in symrhg.F90 routine.
1797 !          XG 001108 :
1798 !          Note that there is a tolerance on the
1799 !          accuracy of tnons, especially when they are found from
1800 !          the symmetry finder (with xred that might be a bit inaccurate)
1801            if (irep > 1) then
1802              do jj=1,nsym_used,irep
1803                setzer=0
1804                ph1r=work2(1,jj);ph1i=work2(2,jj)
1805                do j1=jj,jj+irep-1
1806                  ph2r=work2(1,j1);ph2i=work2(2,j1)
1807                  if (((ph2r+ph1r)**2+(ph2i+ph1i)**2) <= tol14) then
1808                    if (setzer/=1) setzer=-1
1809                  else if (((ph2r-ph1r)**2+(ph2i-ph1i)**2) > tol14) then
1810                    setzer=1
1811                  end if
1812                end do
1813 !              Setzer= 0: phnons are all equal
1814 !              Setzer=-1: phnons are equal in absolute value
1815 !              Setzer= 1: some phnons are different
1816                if (setzer/=0) then
1817                  if (setzer==-1) then
1818                    if (afm_noncoll.and.nspden==4) then
1819                      arg=symafm_used(iperm(jj))
1820                      if (all(symafm_used(iperm(jj:jj+irep-1))==arg)) then
1821                        setzer=1
1822                      else
1823                        do j1=jj,jj+irep-1
1824                          work2(:,j1)=work2(:,j1)*dble(symafm_used(iperm(j1)))
1825                        end do
1826                      end if
1827                    else
1828                      setzer=1
1829                    end if
1830                  end if
1831                  if (setzer==1) work2(:,jj:jj+irep-1)=zero
1832                end if
1833              end do
1834 !            Compress data if irep>1:
1835              jj=0
1836              do isym=1,nsym_used,irep
1837                jj=jj+1
1838                class(jj)=class(isym)
1839                work2(1,jj)=work2(1,isym)
1840                work2(2,jj)=work2(2,isym)
1841              end do
1842            end if
1843 
1844 !          Put new unique points into irrzon array:
1845            irrzon(1+npt:ipt+npt,1,imagn)=class(1:ipt)
1846 
1847 !          Put repetition number into irrzon array:
1848            irrzon(1+nzone,2,imagn)=irep
1849 
1850 !          DEBUG
1851 !          write(std_out,'(a,6i7)' )' irrzg : izone,i1,i2,i3,imagn,irrzon(859,2,1)=',&
1852 !          &      1+nzone,i1,i2,i3,imagn,irrzon(859,2,1)
1853 !          ENDDEBUG
1854 
1855 !          Put phases (or 0) in phnons array:
1856            phnons(:,1+npt:ipt+npt,imagn)=work2(:,1:ipt)
1857 
1858 !          Update number of points in irrzon array:
1859 !          (irep must divide evenly into nsym_used !)
1860            npt=npt+ipt
1861 
1862 !          Update number of classes:
1863            nzone=nzone+1
1864 
1865          end if
1866 !
1867 !        End of loop on reciprocal space points, with indices i1, i2, i3
1868        end do
1869      end do
1870    end do
1871 
1872    if (allocated(symafm_used))  then
1873      ABI_DEALLOCATE(symafm_used)
1874    end if
1875    if (allocated(symrel_used))  then
1876      ABI_DEALLOCATE(symrel_used)
1877    end if
1878    if (allocated(tnons_used))  then
1879      ABI_DEALLOCATE(tnons_used)
1880    end if
1881 
1882  end do ! imagn
1883 
1884 !Make sure number of real space points accounted for equals actual number of grid points
1885  if (npt/=n1*n2*n3) then
1886    write(message, '(a,a,a,a,i10,a,i10,a,a,a,a,a,a,a,a,a)' ) ch10,&
1887 &   ' irrzg : ERROR -',ch10,&
1888 &   '  npt=',npt,' and n1*n2*n3=',n1*n2*n3,' are not equal',ch10,&
1889 &   '  This says that the total of all points in the irreducible',&
1890 &   '  sector in real space',ch10,&
1891 &   '  and all symmetrically equivalent',&
1892 &   '  points, npt, does not equal the actual number',ch10,&
1893 &   '  of real space grid points.'
1894    call wrtout(std_out,message,'COLL')
1895    write(message,'(3a)') &
1896 &   ' This may mean that the input symmetries do not form a group',ch10,&
1897 &   ' Action : check input symmetries carefully for errors.'
1898    MSG_ERROR(message)
1899  end if
1900 
1901 !Perform some checks
1902  do imagn=1,sppoldbl
1903 
1904    do ifft=1,nfftot
1905      if (irrzon(ifft,1,imagn)<1.or.irrzon(ifft,1,imagn)>nfftot) then
1906        write(message,'(a,4i0,a,a)')&
1907 &       '  ifft,irrzon(ifft,1,imagn),nfftot,imagn=',ifft,irrzon(ifft,1,imagn),nfftot,imagn,ch10,&
1908 &       '  =>irrzon goes outside acceptable bounds.'
1909        MSG_BUG(message)
1910      end if
1911    end do
1912 
1913    izonemax=0
1914    do izone=1,nfftot
1915 !    Global bounds
1916      if (irrzon(izone,2,imagn)<0.or.irrzon(izone,2,imagn)>(nsym/imagn)) then
1917        write(message, '(a,5i7,a,a)' )&
1918 &       ' izone,nzone,irrzon(izone,2,imagn),nsym,imagn =',izone,nzone,irrzon(izone,2,imagn),nsym,imagn,ch10,&
1919 &       '  =>irrzon goes outside acceptable bounds.'
1920        MSG_BUG(message)
1921      end if
1922 !    Second index only goes up to nzone
1923      if(izonemax==0)then
1924        if (irrzon(izone,2,imagn)==0)izonemax=izone-1
1925      end if
1926      if(izonemax/=0)then
1927        if (irrzon(izone,2,imagn)/=0) then
1928          message = ' beyond izonemax, irrzon(izone,2,imagn) should be zero'
1929          MSG_BUG(message)
1930        end if
1931      end if
1932    end do
1933 
1934  end do ! imagn
1935 
1936  ABI_DEALLOCATE(class)
1937  ABI_DEALLOCATE(iperm)
1938  ABI_DEALLOCATE(work1)
1939  ABI_DEALLOCATE(work2)
1940 
1941 end subroutine irrzg

m_spacepar/laplacian [ Functions ]

[ Top ] [ m_spacepar ] [ Functions ]

NAME

 laplacian

FUNCTION

 compute the laplacian of a function defined in real space
 the code is written in the way of /3xc/xcden.F90

INPUTS

  gprimd(3,3)=dimensional reciprocal space primitive translations
  mpi_enreg=information about MPI parallelization
  nfft=number of points of the fft grid
  nfunc=number of functions on the grid for which the laplacian is to be calculated
  ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft
  paral_kgb=flag controlling (k,g,bands) parallelization
  (optional) rdfuncr(nfft,nfunc)=real(dp) discretized functions in real space
  rdfuncg_in TO BE DESCRIBED SB 090901
  laplacerdfuncg_in TO BE DESCRIBED SB 090901
  (optional) g2cart_in(nfft) = G**2 on the grid

OUTPUT

 (optional) laplacerdfuncr = laplacian in real space of the functions in rdfuncr
 (optional) rdfuncg = real(dp) discretized functions in fourier space
 (optional) laplacerdfuncg = real(dp) discretized laplacian of the functions in fourier space
 (optional) g2cart_out(nfft) = G**2 on the grid
  rdfuncg_out TO BE DESCRIBED SB 090901
  laplacerdfuncg_out TO BE DESCRIBED SB 090901

PARENTS

      frskerker1,frskerker2,moddiel_csrb,prcrskerker1,prcrskerker2

CHILDREN

      fourdp,ptabs_fourdp

SOURCE

565 subroutine laplacian(gprimd,mpi_enreg,nfft,nfunc,ngfft,paral_kgb,rdfuncr,&
566 &  laplacerdfuncr,rdfuncg_out,laplacerdfuncg_out,g2cart_out,rdfuncg_in,g2cart_in)
567 
568 
569 !This section has been created automatically by the script Abilint (TD).
570 !Do not modify the following lines by hand.
571 #undef ABI_FUNC
572 #define ABI_FUNC 'laplacian'
573 !End of the abilint section
574 
575  implicit none
576 
577 !Arguments ------------------------------------
578 !scalars
579  integer,intent(in) :: nfft,nfunc,paral_kgb
580  type(MPI_type),intent(in) :: mpi_enreg
581 !arrays
582  integer,intent(in) :: ngfft(18)
583  real(dp),intent(in) :: gprimd(3,3)
584  real(dp),intent(inout),optional :: laplacerdfuncr(nfft,nfunc)
585  real(dp),intent(inout),optional,target :: rdfuncr(nfft,nfunc)
586  real(dp),intent(in),optional,target :: g2cart_in(nfft) !vz_i
587  real(dp),intent(out),optional,target :: g2cart_out(nfft)  !vz_i
588  real(dp),intent(out),optional,target :: laplacerdfuncg_out(2,nfft,nfunc)
589  real(dp),intent(in),optional,target :: rdfuncg_in(2,nfft,nfunc) !vz_i
590  real(dp),intent(out),optional,target :: rdfuncg_out(2,nfft,nfunc)
591 
592 !Local variables-------------------------------
593 !scalars
594  integer :: count,i1,i2,i3,id1,id2,id3,ifft,ifunc,ig1,ig2,ig3,ii1,n1,n2
595  integer :: n3
596  real(dp) :: b11,b12,b13,b21,b22,b23,b31,b32,b33
597 !arrays
598  integer, ABI_CONTIGUOUS pointer :: fftn2_distrib(:),ffti2_local(:)
599  integer, ABI_CONTIGUOUS pointer :: fftn3_distrib(:),ffti3_local(:)
600  real(dp),ABI_CONTIGUOUS pointer :: g2cart(:),laplacerdfuncg(:,:,:),rdfuncg(:,:,:)
601 
602 ! *************************************************************************
603 
604 !Keep local copy of fft dimensions
605  n1=ngfft(1)
606  n2=ngfft(2)
607  n3=ngfft(3)
608 
609  if(present(laplacerdfuncg_out)) then
610    laplacerdfuncg => laplacerdfuncg_out
611  else
612    ABI_ALLOCATE(laplacerdfuncg,(2,nfft,nfunc))
613  end if
614 
615  ! Get the distrib associated with this fft_grid
616  call ptabs_fourdp(mpi_enreg,n2,n3,fftn2_distrib,ffti2_local,fftn3_distrib,ffti3_local)
617 
618 !change the real density rdfuncr on real space on the real density
619 !rdfuncg in reciprocal space
620  if(.not.present(rdfuncg_in)) then
621    if(present(rdfuncg_out)) then
622      rdfuncg => rdfuncg_out
623    else
624      ABI_ALLOCATE(rdfuncg,(2,nfft,nfunc))
625    end if
626    if(present(rdfuncr)) then
627      do ifunc=1,nfunc
628        call fourdp(1,rdfuncg(:,:,ifunc),rdfuncr(:,ifunc),-1,mpi_enreg,nfft,ngfft,paral_kgb,0)
629      end do
630    end if
631  else
632    rdfuncg => rdfuncg_in
633  end if
634 
635 !apply the laplacian on laplacerdfuncr
636 !code from /3xc/xcden.F90
637 !see also src/5common/hatre.F90 and src/5common/moddiel.F90
638 !Keep local copy of fft dimensions
639 !Initialize computation of G^2 in cartesian coordinates
640  if(.not.present(g2cart_in)) then
641    if(present(g2cart_out)) then
642      g2cart => g2cart_out
643    else
644      ABI_ALLOCATE(g2cart,(nfft))
645    end if
646    id1=int(n1/2)+2
647    id2=int(n2/2)+2
648    id3=int(n3/2)+2
649    count=0
650    do i3=1,n3
651      ifft=(i3-1)*n1*(n2/mpi_enreg%nproc_fft)
652      ig3=i3-int(i3/id3)*n3-1
653      do i2=1,n2
654        if (fftn2_distrib(i2)==mpi_enreg%me_fft) then
655          ig2=i2-int(i2/id2)*n2-1
656 
657          ii1=1
658          do i1=ii1,n1
659            ig1=i1-int(i1/id1)*n1-1
660            ifft=ifft+1
661 
662            b11=gprimd(1,1)*real(ig1,dp)
663            b21=gprimd(2,1)*real(ig1,dp)
664            b31=gprimd(3,1)*real(ig1,dp)
665            b12=gprimd(1,2)*real(ig2,dp)
666            b22=gprimd(2,2)*real(ig2,dp)
667            b32=gprimd(3,2)*real(ig2,dp)
668            b13=gprimd(1,3)*real(ig3,dp)
669            b23=gprimd(2,3)*real(ig3,dp)
670            b33=gprimd(3,3)*real(ig3,dp)
671 
672            g2cart(ifft)=( &
673 &           (b11+b12+b13)**2&
674 &           +(b21+b22+b23)**2&
675 &           +(b31+b32+b33)**2&
676 &           )
677            do ifunc=1,nfunc
678 !            compute the laplacian in Fourier space that is * (i x 2pi x G)**2
679              laplacerdfuncg(1,ifft,ifunc) = -rdfuncg(1,ifft,ifunc)*g2cart(ifft)*two_pi*two_pi
680              laplacerdfuncg(2,ifft,ifunc) = -rdfuncg(2,ifft,ifunc)*g2cart(ifft)*two_pi*two_pi
681            end do
682          end do
683        end if
684      end do
685    end do
686    if(.not.present(g2cart_out))  then
687      ABI_DEALLOCATE(g2cart)
688    end if
689  else
690    g2cart => g2cart_in
691    do ifunc=1,nfunc
692      do ifft=1,nfft
693 !      compute the laplacian in Fourier space that is * (i x 2pi x G)**2
694        laplacerdfuncg(1,ifft,ifunc) = -rdfuncg(1,ifft,ifunc)*g2cart(ifft)*two_pi*two_pi
695        laplacerdfuncg(2,ifft,ifunc) = -rdfuncg(2,ifft,ifunc)*g2cart(ifft)*two_pi*two_pi
696      end do
697    end do
698  end if
699 
700 !get the result back into real space
701  if(present(laplacerdfuncr)) then
702    do ifunc=1,nfunc
703      call fourdp(1,laplacerdfuncg(:,:,ifunc),laplacerdfuncr(:,ifunc),1,mpi_enreg,nfft,ngfft,paral_kgb,0)
704    end do
705  end if
706 
707 !deallocate pointers
708  if((.not.present(rdfuncg_in)).and.(.not.present(rdfuncg_in)))  then
709    ABI_DEALLOCATE(rdfuncg)
710  end if
711  if(.not.present(laplacerdfuncg_out))  then
712    ABI_DEALLOCATE(laplacerdfuncg)
713  end if
714 
715 end subroutine laplacian

m_spacepar/meanvalue_g [ Functions ]

[ Top ] [ m_spacepar ] [ Functions ]

NAME

 meanvalue_g

FUNCTION

  Compute the mean value of one wavefunction, in reciprocal space,
  for an operator that is real, diagonal in G-space: <wf|op|wf>
  For the time being, only spin-independent operators are treated.

INPUTS

  diag(npw)=diagonal operator (real, spin-independent!)
  filter= if 1, need to filter on the value of diag, that must be less than huge(0.0d0)*1.d-11
          otherwise, should be 0
  istwf_k=storage mode of the vectors
  npw=number of planewaves of the vector
  nspinor=number of spinor components
  vect(2,npw*nspinor)=vector
  vect1(2,npw*nspinor*use_ndo)=vector1 (=vector in most of the cases)
  use_ndo = says if vect=/vect1

OUTPUT

  ar=mean value

PARENTS

      dfpt_vtowfk,energy,forstrnps,vtowfk

CHILDREN

      xmpi_sum

SOURCE

344 subroutine meanvalue_g(ar,diag,filter,istwf_k,mpi_enreg,npw,nspinor,vect,vect1,use_ndo,ar_im)
345 
346 
347 !This section has been created automatically by the script Abilint (TD).
348 !Do not modify the following lines by hand.
349 #undef ABI_FUNC
350 #define ABI_FUNC 'meanvalue_g'
351 !End of the abilint section
352 
353  implicit none
354 
355 !Arguments ------------------------------------
356 !scalars
357  integer,intent(in) :: filter,istwf_k,npw,nspinor,use_ndo
358  real(dp),intent(out) :: ar
359  real(dp),intent(out),optional :: ar_im
360  type(MPI_type),intent(in) :: mpi_enreg
361 !arrays
362  real(dp),intent(in) :: diag(npw),vect(2,npw*nspinor)
363  real(dp),intent(in) :: vect1(2,npw*nspinor)
364 
365 !Local variables-------------------------------
366 !scalars
367  integer :: i1,ierr,ipw,jpw,me_g0
368  character(len=500) :: message
369 !arrays
370 
371 ! *************************************************************************
372  me_g0 = mpi_enreg%me_g0
373 
374  DBG_CHECK(ANY(filter==(/0,1/)),"Wrong filter")
375  DBG_CHECK(ANY(nspinor==(/1,2/)),"Wrong nspinor")
376  DBG_CHECK(ANY(istwf_k==(/(ipw,ipw=1,9)/)),"Wrong istwf_k")
377 
378  if(nspinor==2 .and. istwf_k/=1)then
379    write(message,'(a,a,a,i6,a,i6)')&
380 &   'When istwf_k/=1, nspinor must be 1,',ch10,&
381 &   'however, nspinor=',nspinor,', and istwf_k=',istwf_k
382    MSG_BUG(message)
383  end if
384 
385  if(use_ndo==1 .and. (istwf_k==2 .and.me_g0==1)) then
386    MSG_BUG('use_ndo==1, not tested')
387  end if
388 
389  ar=zero
390  if(present(ar_im)) ar_im=zero
391 
392 !Normal storage mode
393  if(istwf_k==1)then
394 
395 !  No filter
396    if(filter==0)then
397 !$OMP PARALLEL DO REDUCTION(+:ar)
398      do ipw=1,npw
399        ar=ar+diag(ipw)*(vect(1,ipw)*vect1(1,ipw)+vect(2,ipw)*vect1(2,ipw))
400      end do
401      if(nspinor==2)then
402 !$OMP PARALLEL DO REDUCTION(+:ar) PRIVATE(jpw)
403        do ipw=1+npw,2*npw
404          jpw=ipw-npw
405          ar=ar+diag(jpw)*(vect(1,ipw)*vect1(1,ipw)+vect(2,ipw)*vect1(2,ipw))
406        end do
407      end if
408      if(use_ndo==1 .and. nspinor==2)then
409 !$OMP PARALLEL DO REDUCTION(+:ar_im)
410        do ipw=1,npw
411          ar_im=ar_im+diag(ipw)*(vect1(1,ipw)*vect(2,ipw)-vect1(2,ipw)*vect(1,ipw))
412        end do
413 !$OMP PARALLEL DO REDUCTION(+:ar_im) PRIVATE(jpw)
414        do ipw=1+npw,2*npw
415          jpw=ipw-npw
416          ar_im=ar_im+diag(jpw)*(vect1(1,ipw)*vect(2,ipw)-vect1(2,ipw)*vect(1,ipw))
417        end do
418      end if
419 
420 !    !$OMP PARALLEL DO REDUCTION(+:ar,ar_im)
421 !    do ipw=1,npw
422 !    ar=ar+diag(ipw)*(vect(1,ipw)*vect1(1,ipw)+vect(2,ipw)*vect1(2,ipw))
423 !    if(use_ndo==1.and.nspinor==2) ar_im=ar_im+diag(ipw)*(vect1(1,ipw)*vect(2,ipw)-vect1(2,ipw)*vect(1,ipw))
424 !    end do
425 !    if(nspinor==2)then
426 !    !$OMP PARALLEL DO PRIVATE(ipw) REDUCTION(+:ar,ar_im)
427 !    do ipw=1+npw,2*npw
428 !    ar=ar+diag(ipw-npw)*(vect(1,ipw)*vect1(1,ipw)+vect(2,ipw)*vect1(2,ipw))
429 !    if(use_ndo==1.and.nspinor==2) ar_im=ar_im+diag(ipw-npw)*(vect1(1,ipw)*vect(2,ipw)-vect1(2,ipw)*vect(1,ipw))
430 !    end do
431 !    end if
432    else ! will filter
433 
434 !$OMP PARALLEL DO REDUCTION(+:ar)
435      do ipw=1,npw
436        if(diag(ipw)<huge(0.0d0)*1.d-11)then
437          ar=ar+diag(ipw)*(vect(1,ipw)*vect1(1,ipw)+vect(2,ipw)*vect1(2,ipw))
438        end if
439      end do
440      if(nspinor==2)then
441 !$OMP PARALLEL DO REDUCTION(+:ar) PRIVATE(jpw)
442        do ipw=1+npw,2*npw
443          jpw=ipw-npw
444          if(diag(jpw)<huge(0.0d0)*1.d-11)then
445            ar=ar+diag(jpw)*(vect(1,ipw)*vect1(1,ipw)+vect(2,ipw)*vect1(2,ipw))
446          end if
447        end do
448      end if
449      if(use_ndo==1 .and. nspinor==2)then
450 !$OMP PARALLEL DO REDUCTION(+:ar_im)
451        do ipw=1,npw
452          if(diag(ipw)<huge(0.0d0)*1.d-11)then
453            ar_im=ar_im+diag(ipw)*(vect1(1,ipw)*vect(2,ipw)-vect1(2,ipw)*vect(1,ipw))
454          end if
455        end do
456 !$OMP PARALLEL DO REDUCTION(+:ar_im) PRIVATE(jpw)
457        do ipw=1+npw,2*npw
458          jpw=ipw-npw
459          if(diag(jpw)<huge(0.0d0)*1.d-11)then
460            ar_im=ar_im+diag(jpw)*(vect1(1,ipw)*vect(2,ipw)-vect1(2,ipw)*vect(1,ipw))
461          end if
462        end do
463      end if
464 
465 
466 !    !$OMP PARALLEL DO PRIVATE(ipw) REDUCTION(+:ar,ar_im)
467 !    do ipw=1,npw
468 !    if(diag(ipw)<huge(0.0d0)*1.d-11)then
469 !    ar=ar+diag(ipw)*(vect(1,ipw)*vect1(1,ipw)+vect(2,ipw)*vect1(2,ipw))
470 !    if(use_ndo==1.and.nspinor==2) ar_im=ar_im+diag(ipw)*(vect1(1,ipw)*vect(2,ipw)-vect1(2,ipw)*vect(1,ipw))
471 !    end if
472 !    end do
473 !    if(nspinor==2)then
474 !    !$OMP PARALLEL DO PRIVATE(ipw) REDUCTION(+:ar,ar_im)
475 !    do ipw=1+npw,2*npw
476 !    if(diag(ipw-npw)<huge(0.0d0)*1.d-11)then
477 !    ar=ar+diag(ipw-npw)*(vect(1,ipw)*vect1(1,ipw)+vect(2,ipw)*vect1(2,ipw))
478 !    if(use_ndo==1.and.nspinor==2) ar_im=ar_im+diag(ipw-npw)*(vect1(1,ipw)*vect(2,ipw)-vect1(2,ipw)*vect(1,ipw))
479 !    end if
480 !    end do
481 !    end if ! nspinor==2
482 
483    end if ! filter==0
484 
485  else if(istwf_k>=2)then
486 
487    if(filter==0)then
488      i1=1
489      if(istwf_k==2 .and. me_g0==1)then ! MPIWF need to know which proc has G=0
490        ar=half*diag(1)*vect(1,1)*vect1(1,1) ; i1=2
491      end if
492 
493 !$OMP PARALLEL DO REDUCTION(+:ar)
494      do ipw=i1,npw
495        ar=ar+diag(ipw)*(vect(1,ipw)*vect1(1,ipw)+vect(2,ipw)*vect1(2,ipw))
496      end do
497 
498    else ! filter/=0
499      i1=1
500      if(istwf_k==2 .and. me_g0==1)then
501        if(diag(1)<huge(0.0d0)*1.d-11)then
502          ar=half*diag(1)*vect(1,1)*vect1(1,1) ; i1=2
503        end if
504      end if
505 
506 !$OMP PARALLEL DO REDUCTION(+:ar)
507      do ipw=i1,npw
508        if(diag(ipw)<huge(0.0d0)*1.d-11)then
509          ar=ar+diag(ipw)*(vect(1,ipw)*vect1(1,ipw)+vect(2,ipw)*vect1(2,ipw))
510        end if
511      end do
512    end if ! filter==0
513 
514    ar=two*ar
515 
516  end if ! istwf_k
517 
518 !MPIWF need to make reduction on ar and ai .
519  if(mpi_enreg%paral_kgb==1)then
520    call xmpi_sum(ar,mpi_enreg%comm_bandspinorfft ,ierr)
521    if(present(ar_im))then
522      call xmpi_sum(ar_im,mpi_enreg%comm_bandspinorfft,ierr)
523    end if
524  end if
525 
526 end subroutine meanvalue_g

m_spacepar/redgr [ Functions ]

[ Top ] [ m_spacepar ] [ Functions ]

NAME

 redgr

FUNCTION

 Compute reduced gradients of a real function on the usual unshifted
 fft grid. The gradient directions are the along the primitive
 reciprocal lattice vectors.
 The input function is intended to be a single spin component of
 the valence charge density, the valence + core charge densities
 or the first-order core charge density for use in frozen wf
 elastic tensor calculations within the GGA.

NOTES

 Closely linked to xcden, but limited to Q=0, real charge densities,
 and unshifted grids.

INPUTS

  mpi_enreg=information about MPI parallelization
  nfft=(effective) number of FFT grid points (for this processor)
  ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft
  frin(nfft)=real space input function

OUTPUT

  frredgr(nfft,3)= reduced gradient of input function (same units as frin)

PARENTS

      dfpt_eltfrxc

CHILDREN

      fourdp

SOURCE

752 subroutine redgr (frin,frredgr,mpi_enreg,nfft,ngfft,paral_kgb)
753 
754 
755 !This section has been created automatically by the script Abilint (TD).
756 !Do not modify the following lines by hand.
757 #undef ABI_FUNC
758 #define ABI_FUNC 'redgr'
759 !End of the abilint section
760 
761  implicit none
762 
763 !Arguments ------------------------------------
764 !scalars
765  integer,intent(in) :: nfft,paral_kgb
766  type(MPI_type),intent(in) :: mpi_enreg
767 !arrays
768  integer,intent(in) :: ngfft(18)
769  real(dp),intent(in) :: frin(nfft)
770  real(dp),intent(out) :: frredgr(nfft,3)
771 
772 !Local variables-------------------------------
773 !scalars
774  integer :: cplex_tmp,i1,i2,i3,id,idir,ifft,ig,ii,ing,n1,n2,n3
775 !arrays
776  real(dp),allocatable :: gg(:,:),wkcmpx(:,:),work(:),workgr(:,:)
777 
778 ! *************************************************************************
779 
780 !Only real arrays are treated
781  cplex_tmp=1
782 
783 !Keep local copy of fft dimensions
784  n1=ngfft(1) ; n2=ngfft(2) ; n3=ngfft(3)
785 
786 !In order to speed the routine, precompute the components of g, including 2pi factor
787  ABI_ALLOCATE(gg,(max(n1,n2,n3),3))
788  do ii=1,3
789    id=ngfft(ii)/2+2
790    do ing=1,ngfft(ii)
791      ig=ing-(ing/id)*ngfft(ii)-1
792      gg(ing,ii)=two_pi*ig
793    end do
794 !  Note that the G <-> -G symmetry must be maintained
795    if(mod(ngfft(ii),2)==0)gg(ngfft(ii)/2+1,ii)=zero
796  end do
797 
798  ABI_ALLOCATE(wkcmpx,(2,nfft))
799  ABI_ALLOCATE(work,(nfft))
800  ABI_ALLOCATE(workgr,(2,nfft))
801 
802 !Obtain rho(G) in wkcmpx from input rho(r)
803  work(:)=frin(:)
804 
805  call fourdp(cplex_tmp,wkcmpx,work,-1,mpi_enreg,nfft,ngfft,paral_kgb,0)
806 
807 !Gradient calculation for three reduced components in turn.
808 !Code duplicated to remove logic from loops.
809  do idir=1,3
810    if(idir==1) then
811 !$OMP PARALLEL DO PRIVATE(ifft)
812      do i3=1,n3
813        ifft=(i3-1)*n1*n2
814        do i2=1,n2
815          do i1=1,n1
816            ifft=ifft+1
817 !          Multiply by i 2pi G(idir)
818            workgr(2,ifft)= gg(i1,idir)*wkcmpx(1,ifft)
819            workgr(1,ifft)=-gg(i1,idir)*wkcmpx(2,ifft)
820          end do
821        end do
822      end do
823    else if(idir==2) then
824 !$OMP PARALLEL DO PRIVATE(ifft)
825      do i3=1,n3
826        ifft=(i3-1)*n1*n2
827        do i2=1,n2
828          do i1=1,n1
829            ifft=ifft+1
830 !          Multiply by i 2pi G(idir)
831            workgr(2,ifft)= gg(i2,idir)*wkcmpx(1,ifft)
832            workgr(1,ifft)=-gg(i2,idir)*wkcmpx(2,ifft)
833          end do
834        end do
835      end do
836    else
837 !$OMP PARALLEL DO PRIVATE(ifft)
838      do i3=1,n3
839        ifft=(i3-1)*n1*n2
840        do i2=1,n2
841          do i1=1,n1
842            ifft=ifft+1
843 !          Multiply by i 2pi G(idir)
844            workgr(2,ifft)= gg(i3,idir)*wkcmpx(1,ifft)
845            workgr(1,ifft)=-gg(i3,idir)*wkcmpx(2,ifft)
846          end do
847        end do
848      end do
849    end if !idir
850 
851    call fourdp(cplex_tmp,workgr,work,1,mpi_enreg,nfft,ngfft,paral_kgb,0)
852 
853 !$OMP PARALLEL DO
854    do ifft=1,nfft
855      frredgr(ifft,idir)=work(ifft)
856    end do
857 
858  end do !idir
859 
860  ABI_DEALLOCATE(gg)
861  ABI_DEALLOCATE(wkcmpx)
862  ABI_DEALLOCATE(work)
863  ABI_DEALLOCATE(workgr)
864 
865 end subroutine redgr

m_spacepar/rotate_rho [ Functions ]

[ Top ] [ m_spacepar ] [ Functions ]

NAME

 rotate_rho

FUNCTION

 rotate density in real and reciprocal space

INPUTS

  cplex: if 1, real space 1-order functions on FFT grid are REAL, if 2, COMPLEX
  mpi_enreg=informations about MPI parallelization
  nfft=(effective) number of FFT grid points (for this proc) for the "fine" grid (see NOTES in respfn.F90)
  ngfft=array of dimensions for different FFT grids
  nspden=number of spin-density components
  rhor1(cplex*nfft,nspden)=array for Fourier transform of RF electron density
  === if psps%usepaw==1 TODO: extend to PAW
    pawrhoij1(natom) <type(pawrhoij_type)>= 1st-order paw rhoij occupancies and related data
  symrel1=single symmetry operation in real space to apply to rho
  tnon = eventual translation associated to symrel1

OUTPUT

  rhog1_eq(2,nfft)= symmetric density in reciprocal space for equivalent perturbation
  rhor1_eq(cplex*nfft,nspden) = symmetric density in real space for equivalent perturbation

PARENTS

      dfpt_looppert

CHILDREN

      fourdp

TODO

  This routine is deprecated. Now the symmetrization of the **potentials** is done in the m_dvdb
  Moreover one should take linear combinations of symmetrical components

SOURCE

1979 subroutine rotate_rho(cplex, itirev, mpi_enreg, nfft, ngfft, nspden, &
1980 &   rhor1, rhog1_eq, rhor1_eq, symrel1, tnon)
1981 
1982 
1983 !This section has been created automatically by the script Abilint (TD).
1984 !Do not modify the following lines by hand.
1985 #undef ABI_FUNC
1986 #define ABI_FUNC 'rotate_rho'
1987 !End of the abilint section
1988 
1989  implicit none
1990 
1991 !args
1992  integer,intent(in) :: cplex, nfft, nspden, itirev
1993  integer,intent(in) :: ngfft(18)
1994 
1995  integer, intent(in) :: symrel1(3,3)
1996  real(dp),intent(in) :: tnon(3)
1997  real(dp),intent(inout) :: rhor1(cplex*nfft,nspden)
1998 
1999  real(dp),intent(out) :: rhog1_eq(2,nfft)
2000  real(dp),intent(out) :: rhor1_eq(cplex*nfft,nspden)
2001 
2002  type(MPI_type),intent(in) :: mpi_enreg
2003 
2004 ! local vars
2005  integer :: id1, id2, id3
2006  integer :: n1, n2, n3, nd2
2007  integer :: l1, l2, l3
2008  integer :: i1, i2, i3, ind1, ind2
2009  integer :: j1, j2, j3
2010  integer :: k1, k2, k3
2011  integer :: nproc_fft, ispden, me_fft
2012  real(dp) :: arg
2013  logical :: t_tnon_nonzero
2014 
2015  real(dp) :: phnon1(2)
2016  real(dp), allocatable :: workg(:,:), workg_eq(:,:)
2017  character(len=500) :: message
2018 
2019 ! *************************************************************************
2020 
2021  n1=ngfft(1);n2=ngfft(2);n3=ngfft(3);nproc_fft=ngfft(10);me_fft=ngfft(11);nd2=n2/nproc_fft
2022 
2023  id1=n1/2+2
2024  id2=n2/2+2
2025  id3=n3/2+2
2026 
2027  rhog1_eq = zero
2028  rhor1_eq = zero
2029 
2030  if (itirev == 2) then
2031    write (message,'(3a,9I4,1a)') 'using time reversal. ',ch10,'Symrel1 = ', symrel1, ch10
2032  else
2033    write (message,'(3a,9I4,1a)') 'no time reversal. ',ch10,'Symrel1 = ', symrel1, ch10
2034  end if
2035  !call wrtout(std_out, message, 'COLL')
2036 
2037  t_tnon_nonzero = (any(abs(tnon) > tol12))
2038 
2039 ! eventually, for FFT parallelization
2040 ! call ptabs_fourdp(mpi_enreg,n2,n3,fftn2_distrib,ffti2_local,fftn3_distrib,ffti3_local)
2041 
2042  ABI_ALLOCATE(workg,(2,nfft))
2043  ABI_ALLOCATE(workg_eq,(2,nfft))
2044  do ispden = 1, nspden
2045 
2046 ! fft input rhor1 to reciprocal space: uses work* as a buffer
2047    call fourdp(cplex,workg,rhor1(:,ispden),-1,mpi_enreg,nfft,ngfft,mpi_enreg%paral_kgb,0)
2048 
2049 ! below taken from irrzg and setsym
2050 !  Loop over reciprocal space grid points:
2051 !  loop over local points in workg, and get back transform from rhog1,
2052 !  which is presumed complete on each proc
2053    ind1=0
2054    do i3=1,n3
2055      do i2=1,n2
2056 !       if(fftn2_distrib(i2)/=me_fft)  cycle ! this ind is not to be treated by me_fft
2057        do i1=1,n1
2058 
2059          ind1=ind1+1
2060 !       r2=ffti2_local(i2+1) - 1
2061 !       ind=n1*(nd2*i3+r2)+i1+1 !this is ind in the current proc
2062 
2063 !      Get location of G vector (grid point) centered at 0 0 0
2064          l1=i1-(i1/id1)*n1-1
2065          l2=i2-(i2/id2)*n2-1
2066          l3=i3-(i3/id3)*n3-1
2067 
2068 !      Get rotated G vector Gj for each symmetry element
2069 !      -- here we use the TRANSPOSE of symrel1; assuming symrel1 expresses
2070 !      the rotation in real space, the transpose is then appropriate
2071 !      for G space symmetrization (p. 1172d,e of notes, 2 June 1995).
2072          j1=symrel1(1,1)*l1+symrel1(2,1)*l2+symrel1(3,1)*l3
2073          j2=symrel1(1,2)*l1+symrel1(2,2)*l2+symrel1(3,2)*l3
2074          j3=symrel1(1,3)*l1+symrel1(2,3)*l2+symrel1(3,3)*l3
2075 
2076 !      Map into [0,n-1] and then add 1 for array index in [1,n]
2077          k1=1+mod(n1+mod(j1,n1),n1)
2078          k2=1+mod(n2+mod(j2,n2),n2)
2079          k3=1+mod(n3+mod(j3,n3),n3)
2080 
2081 !      Get linear index of rotated point Gj
2082          ind2=k1+n1*((k2-1)+n2*(k3-1))
2083 !       r2=ffti2_local(j2+1) - 1
2084 !       ind=n1*(nd2*j3+r2)+j1+1 !this is ind may be in another proc!!
2085 
2086          phnon1(1) = one
2087          phnon1(2) = zero
2088          if (t_tnon_nonzero) then
2089 !        compute exp(-2*Pi*I*G dot tau) using original G
2090 ! NB: this phase is same as that in irrzg and phnons1, and corresponds to complex conjugate of phase from G to Gj;
2091 ! we use it immediately below, to go _to_ workg(ind1)
2092 ! TODO : replace this with complex powers of exp(2pi tnon(1)) etc...
2093            arg=two_pi*(dble(l1)*tnon(1)+dble(l2)*tnon(2)+dble(l3)*tnon(3))
2094            phnon1(1) = cos(arg)
2095            phnon1(2) =-sin(arg)
2096          end if
2097 
2098 !      rho(Strans*G)=exp(2*Pi*I*(G) dot tau_S) rho(G)
2099          workg_eq (1, ind1) = phnon1(1) * workg(1, ind2) &
2100 &         - phnon1(2) * workg(2, ind2)
2101          workg_eq (2, ind1) = phnon1(1) * workg(2, ind2) &
2102 &         + phnon1(2) * workg(1, ind2)
2103 
2104        end do
2105      end do
2106    end do
2107 
2108 ! accumulate rhog1_eq
2109    if (ispden == 1) rhog1_eq = workg_eq
2110 
2111 ! FFT back to real space to get rhor1_eq
2112 !    Pull out full or spin up density, now symmetrized
2113    call fourdp(cplex,workg_eq,rhor1_eq(:,ispden),1,mpi_enreg,nfft,ngfft,mpi_enreg%paral_kgb,0)
2114 
2115  end do !nspden
2116 
2117  ABI_DEALLOCATE(workg)
2118  ABI_DEALLOCATE(workg_eq)
2119 
2120 end subroutine rotate_rho

m_spacepar/setsym [ Functions ]

[ Top ] [ m_spacepar ] [ Functions ]

NAME

 setsym

FUNCTION

 Set up irreducible zone in  G space by direct calculation.
 Do not call this routine if nsym=1 (only identity symmetry).
 Only indsym and symrec get returned if iscf=0.
 symrec needed to symmetrize coordinate gradients in sygrad.
 (symrec is redundant and could be removed later in favor of symrel)

INPUTS

 iscf=(<= 0  =>non-SCF), >0 => SCF
 natom=number of atoms in unit cell
 nfft=(effective) number of FFT grid points (for this processor)
 ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft
 nspden=number of spin-density components
 nsppol=1 for unpolarized, 2 for spin-polarized
 nsym=number of symmetries in space group (at least 1)
 symafm(nsym)=(anti)ferromagnetic part of symmetry operations
 symrel(3,3,nsym)=symmetry operations in terms of real space primitive translations
 tnons(3,nsym)=nonsymmorphic translations of space group in terms
 of real space primitive translations (may be 0)
 typat(natom)=atom type (integer) for each atom
 xred(3,natom)=atomic coordinates in terms of real space primitive translations

OUTPUT

 indsym(4,nsym,natom)=indirect indexing of atom labels--see subroutine
   symatm for definition (if nsym>1)
 irrzon(nfft,2,(nspden/nsppol)-3*(nspden/4))=irreducible zone data
 phnons(2,nfft,(nspden/nsppol)-3*(nspden/4))=nonsymmorphic translation phases
 symrec(3,3,nsym)=symmetry operations in terms of reciprocal
   space primitive translations (if nsym>1)

NOTES

 nsppol and nspden are needed in case of (anti)ferromagnetic symmetry operations

PARENTS

      dfpt_looppert,gstate,m_dvdb,nonlinear,respfn,scfcv

CHILDREN

      chkgrp,irrzg,mati3inv,symatm,symdet,timab

SOURCE

2168 subroutine setsym(indsym,irrzon,iscf,natom,nfft,ngfft,nspden,nsppol,nsym,phnons,&
2169 & symafm,symrec,symrel,tnons,typat,xred)
2170 
2171 
2172 !This section has been created automatically by the script Abilint (TD).
2173 !Do not modify the following lines by hand.
2174 #undef ABI_FUNC
2175 #define ABI_FUNC 'setsym'
2176 !End of the abilint section
2177 
2178  implicit none
2179 
2180 !Arguments ------------------------------------
2181 !scalars
2182  integer,intent(in) :: iscf,natom,nfft,nspden,nsppol,nsym
2183 !arrays
2184  integer,intent(in) :: ngfft(18),symafm(nsym),symrel(3,3,nsym),typat(natom)
2185  integer,intent(out) :: indsym(4,nsym,natom)
2186  integer,intent(inout) :: irrzon(nfft,2,(nspden/nsppol)-3*(nspden/4)) !vz_i
2187  integer,intent(out) :: symrec(3,3,nsym)
2188  real(dp),intent(in) :: tnons(3,nsym),xred(3,natom)
2189  real(dp),intent(out) :: phnons(2,nfft,(nspden/nsppol)-3*(nspden/4))
2190 
2191 !Local variables-------------------------------
2192 !scalars
2193  integer :: isym,ierr
2194  real(dp) :: tolsym8
2195 !arrays
2196  integer,allocatable :: determinant(:)
2197  real(dp) :: tsec(2)
2198 
2199 ! *************************************************************************
2200 
2201  call timab(6,1,tsec)
2202 
2203 !Check that symmetries have unity determinant
2204  ABI_ALLOCATE(determinant,(nsym))
2205  call symdet(determinant,nsym,symrel)
2206  ABI_DEALLOCATE(determinant)
2207 
2208 
2209 !Get the symmetry matrices in terms of reciprocal basis
2210  do isym=1,nsym
2211    call mati3inv(symrel(:,:,isym),symrec(:,:,isym))
2212  end do
2213 
2214 !Check for group closure
2215  call chkgrp(nsym,symafm,symrel,ierr)
2216  ABI_CHECK(ierr==0,"Error in group closure")
2217 
2218  call chkgrp(nsym,symafm,symrec,ierr)
2219  ABI_CHECK(ierr==0,"Error in group closure")
2220 
2221 !Obtain a list of rotated atom labels:
2222  tolsym8=tol8
2223  call symatm(indsym,natom,nsym,symrec,tnons,tolsym8,typat,xred)
2224 
2225 !If non-SCF calculation, or nsym==1, do not need IBZ data
2226  if ( (iscf>0 .or. iscf==-3) .and. nsym>1 ) then
2227 !  Locate irreducible zone in reciprocal space for symmetrization:
2228    call irrzg(irrzon,nspden,nsppol,nsym,ngfft(1),ngfft(2),ngfft(3),phnons,symafm,symrel,tnons)
2229  end if
2230 
2231  call timab(6,2,tsec)
2232 
2233 end subroutine setsym

m_spacepar/symrhg [ Functions ]

[ Top ] [ m_spacepar ] [ Functions ]

NAME

 symrhg

FUNCTION

 From rho(r), generate rho(G), symmetrize it, and
 come back to the real space for a symmetrized rho(r).

INPUTS

 cplex=1 if rhor is real, 2 if rhor is complex
 gprimd(3,3)=dimensional reciprocal space primitive translations
 irrzon(nfft,2,(nspden/nsppol)-3*(nspden/4))=irreducible zone data
 mpi_enreg=information about MPI parallelization
 nfft=(effective) number of FFT grid points (for this processor)
 ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft
 nspden=number of spin-density components
 nsppol=1 for unpolarized, 2 for spin-polarized
 nsym=number of symmetry elements.
 phnons(2,nfft,(nspden/nsppol)-3*(nspden/4))=nonsymmorphic translation phases
 rprimd(3,3)=dimensional real space primitive translations
 symafm(nsym)=(anti)ferromagnetic part of symmetry operations
 symrel(3,3,nsym)=symmetry matrices in real space (integers)

OUTPUT

 rhog(2,nfft)=symmetrized rho(G) (total) electron density in G space

SIDE EFFECTS

 Input/Output
 rhor(cplex*nfft,nspden)=array for electron density in electrons/bohr**3.
 Input, but also output, if symmetrization is applied.
 Also output if nspden > 1 (change spin components)

NOTES

 When using spin-polarization (nspden==2),
 put total density in first half of rhor array and spin up in second half
 If (nspden=2 and nsppol=2) the density is transformed as  (up,down) => (up+down,up)
 If (nspden=2 and nsppol=1) anti-ferromagnetic symmetry operations
  must be used, such as to transform (2*up) => (up+down,up)
 In spin-polarized, and if there is no symmetry to be
 applied on the system, only the total density is generated in G space

PARENTS

      dfpt_mkrho,dfpt_nstpaw,dfpt_rhofermi,dfpt_vtorho,m_dvdb,mkrho
      suscep_stat,vtorho,vtorhorec,vtorhotf,wfd_mkrho

CHILDREN

      fourdp,matr3inv,ptabs_fourdp,symredcart,timab,xmpi_sum

SOURCE

1116 subroutine symrhg(cplex,gprimd,irrzon,mpi_enreg,nfft,nfftot,ngfft,nspden,nsppol,nsym,paral_kgb,&
1117 &                 phnons,rhog,rhor,rprimd,symafm,symrel)
1118 
1119 
1120 !This section has been created automatically by the script Abilint (TD).
1121 !Do not modify the following lines by hand.
1122 #undef ABI_FUNC
1123 #define ABI_FUNC 'symrhg'
1124 !End of the abilint section
1125 
1126  implicit none
1127 
1128 !Arguments ------------------------------------
1129 !scalars
1130  integer,intent(in) :: cplex,nfft,nfftot,nspden,nsppol,nsym,paral_kgb
1131  type(MPI_type),intent(in) :: mpi_enreg
1132 !arrays
1133  integer,intent(in) :: irrzon(nfftot**(1-1/nsym),2,(nspden/nsppol)-3*(nspden/4)),ngfft(18)
1134  integer,intent(in) :: symafm(nsym),symrel(3,3,nsym)
1135  real(dp),intent(in) :: gprimd(3,3),phnons(2,nfftot**(1-1/nsym),(nspden/nsppol)-3*(nspden/4)),rprimd(3,3)
1136  real(dp),intent(inout) :: rhor(cplex*nfft,nspden)
1137  real(dp),intent(out) :: rhog(2,nfft)
1138 
1139 !Local variables-------------------------------
1140 !scalars
1141  integer :: ier,imagn,ind,ind2,indsy,ispden,isym,iup,izone,izone_max,j,j1,j2,j3,jsym
1142  integer :: k1,k2,k3,l1,l2,l3,me_fft
1143  integer :: n1,n2,n3,nd2,nproc_fft,nspden_eff,nsym_used,numpt,nup
1144  integer :: r2,rep,spaceComm
1145  logical,parameter :: afm_noncoll=.true.  ! TRUE if antiferro symmetries are used in non-collinear magnetism
1146  real(dp) :: magxsu1,magxsu2,magysu1,magysu2,magzsu1,magzsu2,mxi,mxr,myi,myr,mzi,mzr,phi,phr,rhosu1,rhosu2
1147  !character(len=500) :: message
1148 !arrays
1149  integer,allocatable :: isymg(:)
1150  integer, ABI_CONTIGUOUS pointer :: fftn2_distrib(:),ffti2_local(:)
1151  integer, ABI_CONTIGUOUS pointer :: fftn3_distrib(:),ffti3_local(:)
1152  real(dp) :: tsec(2)
1153  real(dp),allocatable :: magngx(:,:),magngy(:,:),magngz(:,:)
1154  real(dp),allocatable :: rhosu1_arr(:),rhosu2_arr(:),work(:)
1155  real(dp),allocatable :: symafm_used(:),symrec_cart(:,:,:),symrel_cart(:,:,:)
1156 
1157 !*************************************************************************
1158 !
1159 !Note the timing channel 17 excludes the different Fourier transforms
1160 
1161  ABI_ALLOCATE(work,(cplex*nfft))
1162 
1163 !Special treatment for spin-polarized case
1164  if(nspden==2 .and. nsppol==2) then
1165 !  When nspden=2 and nsppol=2, put total density in first half
1166 !  of rhor array and spin up in second half  (up,down) => (up+down,up)
1167    call timab(17,1,tsec)
1168    work(:)=rhor(:,1)               ! up => work
1169    rhor(:,1)=rhor(:,1)+rhor(:,2)   ! up+down
1170    rhor(:,2)=work(:)               ! work => up
1171    call timab(17,2,tsec)
1172  end if
1173 
1174 !Special treatment for antiferromagnetism case
1175  if(nspden==2 .and. nsppol==1) then
1176    call timab(17,1,tsec)
1177 !  When nspden=2 and nsppol=1, (2*up) => (2*up,up)
1178 !  Indeed, what was delivered to the present routine is a "total" density,
1179 !  obtained from occupation numbers varying between 0 and 2,
1180 !  but for spin up only potential.
1181    rhor(:,2)=half*rhor(:,1)
1182    call timab(17,2,tsec)
1183  end if
1184 
1185 !Special treatment for non-collinear magnetism case
1186  if(nspden==4) then
1187    call timab(17,1,tsec)
1188 !FR the half factors missing are recovered in dfpt_mkvxc_noncoll and dfpt_accrho
1189    rhor(:,1)=rhor(:,1)+rhor(:,4)     !nup+ndown
1190    rhor(:,2)=rhor(:,2)-rhor(:,1)     !mx (n+mx-n)
1191    rhor(:,3)=rhor(:,3)-rhor(:,1)     !my (n+my-n)
1192    rhor(:,4)=rhor(:,1)-two*rhor(:,4) !mz=n-2ndown
1193    call timab(17,2,tsec)
1194  end if
1195 
1196 
1197  if(nsym==1)then
1198 
1199    if(nspden==2 .and. nsppol==1) then ! There must be at least one anti-ferromagnetic operation
1200      MSG_BUG('In the antiferromagnetic case, nsym cannot be 1')
1201    end if
1202 
1203 !  If not using symmetry, still want total density in G space rho(G).
1204 !  Fourier transform (incl normalization) to get rho(G)
1205    work(:)=rhor(:,1)
1206    call fourdp(cplex,rhog,work,-1,mpi_enreg,nfft,ngfft,paral_kgb,0)
1207  else
1208 
1209 !  Treat either full density, spin-up density or magnetization
1210 !  Note the decrease of ispden to the value 1, in order to finish
1211 !  with rhog of the total density (and not the spin-up density or magnetization)
1212    nspden_eff=nspden;if (nspden==4) nspden_eff=1
1213    do ispden=nspden_eff,1,-1
1214 
1215 !    Prepare the density to be symmetrized, in the reciprocal space
1216      if(nspden==1 .or. nsppol==2 .or. (nspden==4.and.(.not.afm_noncoll)))then
1217        imagn=1
1218        nsym_used=0
1219        do isym=1,nsym
1220          if(symafm(isym)==1)nsym_used=nsym_used+1
1221 !        DEBUG
1222 !        write(std_out,*)' symrhg : isym,symafm(isym)',isym,symafm(isym)
1223 !        ENDDEBUG
1224        end do
1225      else if(nspden==2 .and. nsppol==1)then   ! antiferromagnetic case
1226        imagn=ispden
1227        nsym_used=nsym/ispden
1228      else if (nspden==4) then
1229        imagn=1
1230        nsym_used=nsym/ispden
1231      end if
1232 
1233 !    write(std_out,*)' symrhg : nsym_used=',nsym_used
1234 
1235 !    rhor -fft-> rhog    (rhog is used as work space)
1236 !    Note : it should be possible to reuse rhog in the antiferromagnetic case this would avoid one FFT
1237      work(:)=rhor(:,ispden)
1238      call fourdp(cplex,rhog,work,-1,mpi_enreg,nfft,ngfft,paral_kgb,0)
1239      if (nspden==4) then
1240        ABI_ALLOCATE(magngx,(2,nfft))
1241        ABI_ALLOCATE(magngy,(2,nfft))
1242        ABI_ALLOCATE(magngz,(2,nfft))
1243        work(:)=rhor(:,2)
1244        call fourdp(cplex,magngx,work,-1,mpi_enreg,nfft,ngfft,paral_kgb,0)
1245        work(:)=rhor(:,3)
1246        call fourdp(cplex,magngy,work,-1,mpi_enreg,nfft,ngfft,paral_kgb,0)
1247        work(:)=rhor(:,4)
1248        call fourdp(cplex,magngz,work,-1,mpi_enreg,nfft,ngfft,paral_kgb,0)
1249      end if
1250 
1251 !    Begins the timing here only , to exclude FFTs
1252      call timab(17,1,tsec)
1253 
1254      n1=ngfft(1);n2=ngfft(2);n3=ngfft(3);nproc_fft=ngfft(10);me_fft=ngfft(11);nd2=n2/nproc_fft
1255 
1256 !    Get the distrib associated with this fft_grid
1257      call ptabs_fourdp(mpi_enreg,n2,n3,fftn2_distrib,ffti2_local,fftn3_distrib,ffti3_local)
1258 
1259 !    The following is only valid for total, up or dn density
1260 !    -------------------------------------------------------
1261 
1262 !    Get maxvalue of izone
1263      izone_max=count(irrzon(:,2,imagn)>0)
1264      ABI_ALLOCATE(rhosu1_arr,(izone_max))
1265      ABI_ALLOCATE(rhosu2_arr,(izone_max))
1266 
1267      numpt=0
1268      do izone=1,nfftot
1269 
1270 !      Get repetition number
1271        rep=irrzon(izone,2,imagn)
1272        if(rep==0)exit
1273 
1274 !      Compute number of unique points in this symm class:
1275        nup=nsym_used/rep
1276 
1277 !      Accumulate charge over equivalent points
1278        rhosu1=zero
1279        rhosu2=zero
1280        do iup=1,nup
1281          ind=irrzon(iup+numpt,1,imagn)
1282          j=ind-1;j1=modulo(j,n1);j2=modulo(j/n1,n2);j3=j/(n1*n2);
1283          if(fftn2_distrib(j2+1)==me_fft)  then ! this ind is to be treated by me_fft
1284            r2=ffti2_local(j2+1) - 1
1285            ind=n1*(nd2*j3+r2)+j1+1 !this is ind in the current proc
1286            rhosu1=rhosu1+rhog(1,ind)*phnons(1,iup+numpt,imagn)&
1287 &           -rhog(2,ind)*phnons(2,iup+numpt,imagn)
1288            rhosu2=rhosu2+rhog(2,ind)*phnons(1,iup+numpt,imagn)&
1289 &           +rhog(1,ind)*phnons(2,iup+numpt,imagn)
1290          end if
1291 
1292        end do
1293        rhosu1=rhosu1/dble(nup)
1294        rhosu2=rhosu2/dble(nup)
1295        rhosu1_arr(izone)=rhosu1
1296        rhosu2_arr(izone)=rhosu2
1297 !      Keep index of how many points have been considered:
1298        numpt=numpt+nup
1299 
1300      end do  ! End loop over izone
1301 
1302 !    Reduction in case of FFT parallelization
1303      if(mpi_enreg%nproc_fft>1)then
1304        spaceComm=mpi_enreg%comm_fft
1305        call xmpi_sum(rhosu1_arr,spaceComm,ier)
1306        call xmpi_sum(rhosu2_arr,spaceComm,ier)
1307      end if
1308 
1309 !    Now symmetrize the density
1310      numpt=0
1311      do izone=1,nfftot
1312 
1313 !      Get repetition number
1314        rep=irrzon(izone,2,imagn)
1315        if(rep==0)exit
1316 
1317 !      Compute number of unique points in this symm class:
1318        nup=nsym_used/rep
1319 
1320 !      Define symmetrized rho(G) at equivalent points:
1321        do iup=1,nup
1322          ind=irrzon(iup+numpt,1,imagn)
1323 !        decompose ind-1=n1(n2 j3+ j2)+j1
1324          j=ind-1;j1=modulo(j,n1);j2=modulo(j/n1,n2);j3=j/(n1*n2);
1325          if(fftn2_distrib(j2+1)==me_fft)  then ! this ind is to be treated by me_fft
1326            r2=ffti2_local(j2+1) - 1
1327 !          ind in the proc ind-1=n1(nd2 j3+ r2)+j1
1328            ind=n1*(nd2*j3+r2)+j1+1 !this is ind in the current proc
1329            rhog(1,ind)=rhosu1_arr(izone)*phnons(1,iup+numpt,imagn)&
1330 &           +rhosu2_arr(izone)*phnons(2,iup+numpt,imagn)
1331            rhog(2,ind)=rhosu2_arr(izone)*phnons(1,iup+numpt,imagn)&
1332 &           -rhosu1_arr(izone)*phnons(2,iup+numpt,imagn)
1333          end if
1334        end do
1335 
1336 !      Keep index of how many points have been considered:
1337        numpt=numpt+nup
1338 
1339      end do ! End loop over izone
1340 
1341      ABI_DEALLOCATE(rhosu1_arr)
1342      ABI_DEALLOCATE(rhosu2_arr)
1343 
1344 !    The following is only valid for magnetization
1345 !    ---------------------------------------------
1346      if (nspden==4) then
1347 
1348 !      Transfer symmetries in cartesian coordinates
1349 !      Compute symmetries in reciprocal space in cartesian coordinates
1350        ABI_ALLOCATE(symrec_cart,(3,3,nsym_used))
1351        ABI_ALLOCATE(symrel_cart,(3,3,nsym_used))
1352        ABI_ALLOCATE(symafm_used,(nsym_used))
1353        jsym=0
1354        do isym=1,nsym
1355          if (symafm(isym)/=1.and.(.not.afm_noncoll)) cycle
1356          jsym=jsym+1
1357          symafm_used(jsym)=dble(symafm(isym))
1358          call symredcart(rprimd,gprimd,symrel_cart(:,:,jsym),symrel(:,:,isym))
1359          call matr3inv(symrel_cart(:,:,jsym),symrec_cart(:,:,jsym))
1360        end do
1361 
1362        numpt=count(irrzon(:,1,imagn)>0)
1363        ABI_ALLOCATE(isymg,(numpt))
1364        isymg=0
1365        ABI_ALLOCATE(rhosu1_arr,(3*izone_max))
1366        ABI_ALLOCATE(rhosu2_arr,(3*izone_max))
1367 
1368 !      Accumulate magnetization over equivalent points
1369 !      Use all symmetries (not only those linking different g points)
1370 !      Use Inverse[Transpose[symrel]]=symrec
1371        numpt=0
1372        do izone=1,izone_max
1373          magxsu1=zero;magxsu2=zero
1374          magysu1=zero;magysu2=zero
1375          magzsu1=zero;magzsu2=zero
1376          ind=irrzon(1+numpt,1,1)
1377          rep=irrzon(izone,2,1)
1378          nup=nsym_used/rep
1379          j=ind-1;l1=modulo(j,n1);l2=modulo(j/n1,n2);l3=j/(n1*n2)
1380          jsym=0
1381          do isym=1,nsym
1382            if (symafm(isym)/=1.and.(.not.afm_noncoll)) cycle
1383            jsym=jsym+1
1384            j1=symrel(1,1,isym)*l1+symrel(2,1,isym)*l2+symrel(3,1,isym)*l3
1385            j2=symrel(1,2,isym)*l1+symrel(2,2,isym)*l2+symrel(3,2,isym)*l3
1386            j3=symrel(1,3,isym)*l1+symrel(2,3,isym)*l2+symrel(3,3,isym)*l3
1387            k1=map_symrhg(j1,n1);k2=map_symrhg(j2,n2);k3=map_symrhg(j3,n3)
1388            indsy=1+k1+n1*(k2+n2*k3)
1389            ind2=-1;iup=numpt
1390            do while (ind2/=indsy.and.iup<numpt+nup)
1391              iup=iup+1;ind2=irrzon(iup,1,1)
1392            end do
1393            if (ind2/=indsy) then
1394              MSG_ERROR("ind2/=indsy in symrhg !")
1395            end if
1396            if (isymg(iup)==0) isymg(iup)=jsym
1397            if(fftn2_distrib(modulo((indsy-1)/n1,n2) + 1) == me_fft ) then  ! this is indsy is to be treated by me_fft
1398              indsy=n1*(nd2*k3+ ffti2_local(k2+1) -1)+k1+1        ! this is indsy in the current proc
1399              phr=phnons(1,iup,imagn);if (rep==1) phr=phr*symafm_used(jsym) !if rep==2, symafm is already included in phnons
1400              phi=phnons(2,iup,imagn);if (rep==1) phi=phi*symafm_used(jsym) !(see irrzg.F90)
1401              mxr=symrel_cart(1,1,jsym)*magngx(1,indsy)+symrel_cart(1,2,jsym)*magngy(1,indsy)+symrel_cart(1,3,jsym)*magngz(1,indsy)
1402              mxi=symrel_cart(1,1,jsym)*magngx(2,indsy)+symrel_cart(1,2,jsym)*magngy(2,indsy)+symrel_cart(1,3,jsym)*magngz(2,indsy)
1403              myr=symrel_cart(2,1,jsym)*magngx(1,indsy)+symrel_cart(2,2,jsym)*magngy(1,indsy)+symrel_cart(2,3,jsym)*magngz(1,indsy)
1404              myi=symrel_cart(2,1,jsym)*magngx(2,indsy)+symrel_cart(2,2,jsym)*magngy(2,indsy)+symrel_cart(2,3,jsym)*magngz(2,indsy)
1405              mzr=symrel_cart(3,1,jsym)*magngx(1,indsy)+symrel_cart(3,2,jsym)*magngy(1,indsy)+symrel_cart(3,3,jsym)*magngz(1,indsy)
1406              mzi=symrel_cart(3,1,jsym)*magngx(2,indsy)+symrel_cart(3,2,jsym)*magngy(2,indsy)+symrel_cart(3,3,jsym)*magngz(2,indsy)
1407              magxsu1=magxsu1+mxr*phr-mxi*phi;magxsu2=magxsu2+mxi*phr+mxr*phi
1408              magysu1=magysu1+myr*phr-myi*phi;magysu2=magysu2+myi*phr+myr*phi
1409              magzsu1=magzsu1+mzr*phr-mzi*phi;magzsu2=magzsu2+mzi*phr+mzr*phi
1410            end if
1411          end do
1412          rhosu1_arr(3*izone-2)=magxsu1/dble(nsym_used)
1413          rhosu1_arr(3*izone-1)=magysu1/dble(nsym_used)
1414          rhosu1_arr(3*izone  )=magzsu1/dble(nsym_used)
1415          rhosu2_arr(3*izone-2)=magxsu2/dble(nsym_used)
1416          rhosu2_arr(3*izone-1)=magysu2/dble(nsym_used)
1417          rhosu2_arr(3*izone  )=magzsu2/dble(nsym_used)
1418          numpt=numpt+nup
1419        end do
1420 
1421 !      Reduction in case of FFT parallelization
1422        if(mpi_enreg%nproc_fft>1)then
1423          spaceComm=mpi_enreg%comm_fft
1424          call xmpi_sum(rhosu1_arr,spaceComm,ier)
1425          call xmpi_sum(rhosu2_arr,spaceComm,ier)
1426        end if
1427 
1428 !      Now symmetrize the magnetization at equivalent points
1429 !      Use Transpose[symrel]
1430        numpt=0
1431        do izone=1,izone_max
1432          rep=irrzon(izone,2,imagn)
1433          nup=nsym_used/rep
1434          do iup=1,nup
1435            ind=irrzon(iup+numpt,1,imagn)
1436            j=ind-1;j1=modulo(j,n1);j2=modulo(j/n1,n2);j3=j/(n1*n2)
1437            if(fftn2_distrib(j2+1)==me_fft)  then ! this ind is to be treated by me_fft
1438              r2=ffti2_local(j2+1) - 1
1439              ind=n1*(nd2*j3+r2)+j1+1  ! this is ind in the current proc
1440              jsym=isymg(iup+numpt)
1441              if (jsym==0) then
1442                MSG_ERROR("jsym=0 in symrhg !")
1443              end if
1444              magxsu1=rhosu1_arr(3*izone-2);magxsu2=rhosu2_arr(3*izone-2)
1445              magysu1=rhosu1_arr(3*izone-1);magysu2=rhosu2_arr(3*izone-1)
1446              magzsu1=rhosu1_arr(3*izone  );magzsu2=rhosu2_arr(3*izone  )
1447              phr=phnons(1,iup,imagn);if (rep==1) phr=phr*symafm_used(jsym) !if rep==2, symafm is already included in phnons
1448              phi=phnons(2,iup,imagn);if (rep==1) phi=phi*symafm_used(jsym) !(see irrzg.F90)
1449              mxr=symrec_cart(1,1,jsym)*magxsu1+symrec_cart(2,1,jsym)*magysu1+symrec_cart(3,1,jsym)*magzsu1
1450              mxi=symrec_cart(1,1,jsym)*magxsu2+symrec_cart(2,1,jsym)*magysu2+symrec_cart(3,1,jsym)*magzsu2
1451              myr=symrec_cart(1,2,jsym)*magxsu1+symrec_cart(2,2,jsym)*magysu1+symrec_cart(3,2,jsym)*magzsu1
1452              myi=symrec_cart(1,2,jsym)*magxsu2+symrec_cart(2,2,jsym)*magysu2+symrec_cart(3,2,jsym)*magzsu2
1453              mzr=symrec_cart(1,3,jsym)*magxsu1+symrec_cart(2,3,jsym)*magysu1+symrec_cart(3,3,jsym)*magzsu1
1454              mzi=symrec_cart(1,3,jsym)*magxsu2+symrec_cart(2,3,jsym)*magysu2+symrec_cart(3,3,jsym)*magzsu2
1455              magngx(1,ind)=mxr*phr+mxi*phi
1456              magngx(2,ind)=mxi*phr-mxr*phi
1457              magngy(1,ind)=myr*phr+myi*phi
1458              magngy(2,ind)=myi*phr-myr*phi
1459              magngz(1,ind)=mzr*phr+mzi*phi
1460              magngz(2,ind)=mzi*phr-mzr*phi
1461            end if
1462          end do
1463          numpt=numpt+nup
1464        end do
1465        ABI_DEALLOCATE(isymg)
1466        ABI_DEALLOCATE(rhosu1_arr)
1467        ABI_DEALLOCATE(rhosu2_arr)
1468        ABI_DEALLOCATE(symrec_cart)
1469        ABI_DEALLOCATE(symrel_cart)
1470        ABI_DEALLOCATE(symafm_used)
1471 
1472      end if ! nspden==4
1473 
1474      call timab(17,2,tsec)
1475 
1476 !    Pull out full or spin up density, now symmetrized
1477      call fourdp(cplex,rhog,work,1,mpi_enreg,nfft,ngfft,paral_kgb,0)
1478      rhor(:,ispden)=work(:)
1479      if (nspden==4) then
1480        call fourdp(cplex,magngx,work,1,mpi_enreg,nfft,ngfft,paral_kgb,0)
1481        rhor(:,2)=work(:)
1482        call fourdp(cplex,magngy,work,1,mpi_enreg,nfft,ngfft,paral_kgb,0)
1483        rhor(:,3)=work(:)
1484        call fourdp(cplex,magngz,work,1,mpi_enreg,nfft,ngfft,paral_kgb,0)
1485        rhor(:,4)=work(:)
1486        ABI_DEALLOCATE(magngx)
1487        ABI_DEALLOCATE(magngy)
1488        ABI_DEALLOCATE(magngz)
1489      end if
1490 
1491    end do ! ispden
1492 
1493  end if !  End on the condition nsym==1
1494 
1495  ABI_DEALLOCATE(work)
1496 
1497  contains
1498 
1499    function map_symrhg(j1,n1)
1500 
1501 
1502 !This section has been created automatically by the script Abilint (TD).
1503 !Do not modify the following lines by hand.
1504 #undef ABI_FUNC
1505 #define ABI_FUNC 'map_symrhg'
1506 !End of the abilint section
1507 
1508    integer :: map_symrhg
1509    integer,intent(in) :: j1,n1
1510    map_symrhg=mod(n1+mod(j1,n1),n1)
1511  end function map_symrhg
1512 
1513 end subroutine symrhg