TABLE OF CONTENTS


ABINIT/symrhg [ Functions ]

[ Top ] [ Functions ]

NAME

 symrhg

FUNCTION

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

COPYRIGHT

 Copyright (C) 1998-2018 ABINIT group (DCA, XG, GMR)
 This file is distributed under the terms of the
 GNU General Public License, see ~abinit/COPYING
 or http://www.gnu.org/copyleft/gpl.txt .
 For the initials of contributors, see ~abinit/doc/developers/contributors.txt .

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

 59 #if defined HAVE_CONFIG_H
 60 #include "config.h"
 61 #endif
 62 
 63 #include "abi_common.h"
 64 
 65 subroutine symrhg(cplex,gprimd,irrzon,mpi_enreg,nfft,nfftot,ngfft,nspden,nsppol,nsym,paral_kgb,&
 66 &                 phnons,rhog,rhor,rprimd,symafm,symrel)
 67 
 68  use defs_basis
 69  use defs_abitypes
 70  use m_errors
 71  use m_profiling_abi
 72  use m_xmpi
 73 
 74  use m_mpinfo,   only : ptabs_fourdp
 75 
 76 !This section has been created automatically by the script Abilint (TD).
 77 !Do not modify the following lines by hand.
 78 #undef ABI_FUNC
 79 #define ABI_FUNC 'symrhg'
 80  use interfaces_18_timing
 81  use interfaces_32_util
 82  use interfaces_41_geometry
 83  use interfaces_53_ffts
 84 !End of the abilint section
 85 
 86  implicit none
 87 
 88 !Arguments ------------------------------------
 89 !scalars
 90  integer,intent(in) :: cplex,nfft,nfftot,nspden,nsppol,nsym,paral_kgb
 91  type(MPI_type),intent(in) :: mpi_enreg
 92 !arrays
 93  integer,intent(in) :: irrzon(nfftot**(1-1/nsym),2,(nspden/nsppol)-3*(nspden/4)),ngfft(18)
 94  integer,intent(in) :: symafm(nsym),symrel(3,3,nsym)
 95  real(dp),intent(in) :: gprimd(3,3),phnons(2,nfftot**(1-1/nsym),(nspden/nsppol)-3*(nspden/4)),rprimd(3,3)
 96  real(dp),intent(inout) :: rhor(cplex*nfft,nspden)
 97  real(dp),intent(out) :: rhog(2,nfft)
 98 
 99 !Local variables-------------------------------
100 !scalars
101  integer :: ier,imagn,ind,ind2,indsy,ispden,isym,iup,izone,izone_max,j,j1,j2,j3,jsym
102  integer :: k1,k2,k3,l1,l2,l3,me_fft
103  integer :: n1,n2,n3,nd2,nproc_fft,nspden_eff,nsym_used,numpt,nup
104  integer :: r2,rep,spaceComm
105  logical,parameter :: afm_noncoll=.true.  ! TRUE if antiferro symmetries are used in non-collinear magnetism
106  real(dp) :: magxsu1,magxsu2,magysu1,magysu2,magzsu1,magzsu2,mxi,mxr,myi,myr,mzi,mzr,phi,phr,rhosu1,rhosu2
107  !character(len=500) :: message
108 !arrays
109  integer,allocatable :: isymg(:)
110  integer, ABI_CONTIGUOUS pointer :: fftn2_distrib(:),ffti2_local(:)
111  integer, ABI_CONTIGUOUS pointer :: fftn3_distrib(:),ffti3_local(:)
112  real(dp) :: tsec(2)
113  real(dp),allocatable :: magngx(:,:),magngy(:,:),magngz(:,:)
114  real(dp),allocatable :: rhosu1_arr(:),rhosu2_arr(:),work(:)
115  real(dp),allocatable :: symafm_used(:),symrec_cart(:,:,:),symrel_cart(:,:,:)
116 
117 !*************************************************************************
118 !
119 !Note the timing channel 17 excludes the different Fourier transforms
120 
121  ABI_ALLOCATE(work,(cplex*nfft))
122 
123 !Special treatment for spin-polarized case
124  if(nspden==2 .and. nsppol==2) then
125 !  When nspden=2 and nsppol=2, put total density in first half
126 !  of rhor array and spin up in second half  (up,down) => (up+down,up)
127    call timab(17,1,tsec)
128    work(:)=rhor(:,1)               ! up => work
129    rhor(:,1)=rhor(:,1)+rhor(:,2)   ! up+down
130    rhor(:,2)=work(:)               ! work => up
131    call timab(17,2,tsec)
132  end if
133 
134 !Special treatment for antiferromagnetism case
135  if(nspden==2 .and. nsppol==1) then
136    call timab(17,1,tsec)
137 !  When nspden=2 and nsppol=1, (2*up) => (2*up,up)
138 !  Indeed, what was delivered to the present routine is a "total" density,
139 !  obtained from occupation numbers varying between 0 and 2,
140 !  but for spin up only potential.
141    rhor(:,2)=half*rhor(:,1)
142    call timab(17,2,tsec)
143  end if
144 
145 !Special treatment for non-collinear magnetism case
146  if(nspden==4) then
147    call timab(17,1,tsec)
148 !FR the half factors missing are recovered in dfpt_mkvxc_noncoll and dfpt_accrho
149    rhor(:,1)=rhor(:,1)+rhor(:,4)     !nup+ndown
150    rhor(:,2)=rhor(:,2)-rhor(:,1)     !mx (n+mx-n)
151    rhor(:,3)=rhor(:,3)-rhor(:,1)     !my (n+my-n)
152    rhor(:,4)=rhor(:,1)-two*rhor(:,4) !mz=n-2ndown
153    call timab(17,2,tsec)
154  end if
155 
156 
157  if(nsym==1)then
158 
159    if(nspden==2 .and. nsppol==1) then ! There must be at least one anti-ferromagnetic operation
160      MSG_BUG('In the antiferromagnetic case, nsym cannot be 1')
161    end if
162 
163 !  If not using symmetry, still want total density in G space rho(G).
164 !  Fourier transform (incl normalization) to get rho(G)
165    work(:)=rhor(:,1)
166    call fourdp(cplex,rhog,work,-1,mpi_enreg,nfft,ngfft,paral_kgb,0)
167  else
168 
169 !  Treat either full density, spin-up density or magnetization
170 !  Note the decrease of ispden to the value 1, in order to finish
171 !  with rhog of the total density (and not the spin-up density or magnetization)
172    nspden_eff=nspden;if (nspden==4) nspden_eff=1
173    do ispden=nspden_eff,1,-1
174 
175 !    Prepare the density to be symmetrized, in the reciprocal space
176      if(nspden==1 .or. nsppol==2 .or. (nspden==4.and.(.not.afm_noncoll)))then
177        imagn=1
178        nsym_used=0
179        do isym=1,nsym
180          if(symafm(isym)==1)nsym_used=nsym_used+1
181 !        DEBUG
182 !        write(std_out,*)' symrhg : isym,symafm(isym)',isym,symafm(isym)
183 !        ENDDEBUG
184        end do
185      else if(nspden==2 .and. nsppol==1)then   ! antiferromagnetic case
186        imagn=ispden
187        nsym_used=nsym/ispden
188      else if (nspden==4) then
189        imagn=1
190        nsym_used=nsym/ispden
191      end if
192 
193 !    write(std_out,*)' symrhg : nsym_used=',nsym_used
194 
195 !    rhor -fft-> rhog    (rhog is used as work space)
196 !    Note : it should be possible to reuse rhog in the antiferromagnetic case this would avoid one FFT
197      work(:)=rhor(:,ispden)
198      call fourdp(cplex,rhog,work,-1,mpi_enreg,nfft,ngfft,paral_kgb,0)
199      if (nspden==4) then
200        ABI_ALLOCATE(magngx,(2,nfft))
201        ABI_ALLOCATE(magngy,(2,nfft))
202        ABI_ALLOCATE(magngz,(2,nfft))
203        work(:)=rhor(:,2)
204        call fourdp(cplex,magngx,work,-1,mpi_enreg,nfft,ngfft,paral_kgb,0)
205        work(:)=rhor(:,3)
206        call fourdp(cplex,magngy,work,-1,mpi_enreg,nfft,ngfft,paral_kgb,0)
207        work(:)=rhor(:,4)
208        call fourdp(cplex,magngz,work,-1,mpi_enreg,nfft,ngfft,paral_kgb,0)
209      end if
210 
211 !    Begins the timing here only , to exclude FFTs
212      call timab(17,1,tsec)
213 
214      n1=ngfft(1);n2=ngfft(2);n3=ngfft(3);nproc_fft=ngfft(10);me_fft=ngfft(11);nd2=n2/nproc_fft
215 
216 !    Get the distrib associated with this fft_grid
217      call ptabs_fourdp(mpi_enreg,n2,n3,fftn2_distrib,ffti2_local,fftn3_distrib,ffti3_local)
218 
219 !    The following is only valid for total, up or dn density
220 !    -------------------------------------------------------
221 
222 !    Get maxvalue of izone
223      izone_max=count(irrzon(:,2,imagn)>0)
224      ABI_ALLOCATE(rhosu1_arr,(izone_max))
225      ABI_ALLOCATE(rhosu2_arr,(izone_max))
226 
227      numpt=0
228      do izone=1,nfftot
229 
230 !      Get repetition number
231        rep=irrzon(izone,2,imagn)
232        if(rep==0)exit
233 
234 !      Compute number of unique points in this symm class:
235        nup=nsym_used/rep
236 
237 !      Accumulate charge over equivalent points
238        rhosu1=zero
239        rhosu2=zero
240        do iup=1,nup
241          ind=irrzon(iup+numpt,1,imagn)
242          j=ind-1;j1=modulo(j,n1);j2=modulo(j/n1,n2);j3=j/(n1*n2);
243          if(fftn2_distrib(j2+1)==me_fft)  then ! this ind is to be treated by me_fft
244            r2=ffti2_local(j2+1) - 1
245            ind=n1*(nd2*j3+r2)+j1+1 !this is ind in the current proc
246            rhosu1=rhosu1+rhog(1,ind)*phnons(1,iup+numpt,imagn)&
247 &           -rhog(2,ind)*phnons(2,iup+numpt,imagn)
248            rhosu2=rhosu2+rhog(2,ind)*phnons(1,iup+numpt,imagn)&
249 &           +rhog(1,ind)*phnons(2,iup+numpt,imagn)
250          end if
251 
252        end do
253        rhosu1=rhosu1/dble(nup)
254        rhosu2=rhosu2/dble(nup)
255        rhosu1_arr(izone)=rhosu1
256        rhosu2_arr(izone)=rhosu2
257 !      Keep index of how many points have been considered:
258        numpt=numpt+nup
259 
260      end do  ! End loop over izone
261 
262 !    Reduction in case of FFT parallelization
263      if(mpi_enreg%nproc_fft>1)then
264        spaceComm=mpi_enreg%comm_fft
265        call xmpi_sum(rhosu1_arr,spaceComm,ier)
266        call xmpi_sum(rhosu2_arr,spaceComm,ier)
267      end if
268 
269 !    Now symmetrize the density
270      numpt=0
271      do izone=1,nfftot
272 
273 !      Get repetition number
274        rep=irrzon(izone,2,imagn)
275        if(rep==0)exit
276 
277 !      Compute number of unique points in this symm class:
278        nup=nsym_used/rep
279 
280 !      Define symmetrized rho(G) at equivalent points:
281        do iup=1,nup
282          ind=irrzon(iup+numpt,1,imagn)
283 !        decompose ind-1=n1(n2 j3+ j2)+j1
284          j=ind-1;j1=modulo(j,n1);j2=modulo(j/n1,n2);j3=j/(n1*n2);
285          if(fftn2_distrib(j2+1)==me_fft)  then ! this ind is to be treated by me_fft
286            r2=ffti2_local(j2+1) - 1
287 !          ind in the proc ind-1=n1(nd2 j3+ r2)+j1
288            ind=n1*(nd2*j3+r2)+j1+1 !this is ind in the current proc
289            rhog(1,ind)=rhosu1_arr(izone)*phnons(1,iup+numpt,imagn)&
290 &           +rhosu2_arr(izone)*phnons(2,iup+numpt,imagn)
291            rhog(2,ind)=rhosu2_arr(izone)*phnons(1,iup+numpt,imagn)&
292 &           -rhosu1_arr(izone)*phnons(2,iup+numpt,imagn)
293          end if
294        end do
295 
296 !      Keep index of how many points have been considered:
297        numpt=numpt+nup
298 
299      end do ! End loop over izone
300 
301      ABI_DEALLOCATE(rhosu1_arr)
302      ABI_DEALLOCATE(rhosu2_arr)
303 
304 !    The following is only valid for magnetization
305 !    ---------------------------------------------
306      if (nspden==4) then
307 
308 !      Transfer symmetries in cartesian coordinates
309 !      Compute symmetries in reciprocal space in cartesian coordinates
310        ABI_ALLOCATE(symrec_cart,(3,3,nsym_used))
311        ABI_ALLOCATE(symrel_cart,(3,3,nsym_used))
312        ABI_ALLOCATE(symafm_used,(nsym_used))
313        jsym=0
314        do isym=1,nsym
315          if (symafm(isym)/=1.and.(.not.afm_noncoll)) cycle
316          jsym=jsym+1
317          symafm_used(jsym)=dble(symafm(isym))
318          call symredcart(rprimd,gprimd,symrel_cart(:,:,jsym),symrel(:,:,isym))
319          call matr3inv(symrel_cart(:,:,jsym),symrec_cart(:,:,jsym))
320        end do
321 
322        numpt=count(irrzon(:,1,imagn)>0)
323        ABI_ALLOCATE(isymg,(numpt))
324        isymg=0
325        ABI_ALLOCATE(rhosu1_arr,(3*izone_max))
326        ABI_ALLOCATE(rhosu2_arr,(3*izone_max))
327 
328 !      Accumulate magnetization over equivalent points
329 !      Use all symmetries (not only those linking different g points)
330 !      Use Inverse[Transpose[symrel]]=symrec
331        numpt=0
332        do izone=1,izone_max
333          magxsu1=zero;magxsu2=zero
334          magysu1=zero;magysu2=zero
335          magzsu1=zero;magzsu2=zero
336          ind=irrzon(1+numpt,1,1)
337          rep=irrzon(izone,2,1)
338          nup=nsym_used/rep
339          j=ind-1;l1=modulo(j,n1);l2=modulo(j/n1,n2);l3=j/(n1*n2)
340          jsym=0
341          do isym=1,nsym
342            if (symafm(isym)/=1.and.(.not.afm_noncoll)) cycle
343            jsym=jsym+1
344            j1=symrel(1,1,isym)*l1+symrel(2,1,isym)*l2+symrel(3,1,isym)*l3
345            j2=symrel(1,2,isym)*l1+symrel(2,2,isym)*l2+symrel(3,2,isym)*l3
346            j3=symrel(1,3,isym)*l1+symrel(2,3,isym)*l2+symrel(3,3,isym)*l3
347            k1=map_symrhg(j1,n1);k2=map_symrhg(j2,n2);k3=map_symrhg(j3,n3)
348            indsy=1+k1+n1*(k2+n2*k3)
349            ind2=-1;iup=numpt
350            do while (ind2/=indsy.and.iup<numpt+nup)
351              iup=iup+1;ind2=irrzon(iup,1,1)
352            end do
353            if (ind2/=indsy) then
354              MSG_ERROR("ind2/=indsy in symrhg !")
355            end if
356            if (isymg(iup)==0) isymg(iup)=jsym
357            if(fftn2_distrib(modulo((indsy-1)/n1,n2) + 1) == me_fft ) then  ! this is indsy is to be treated by me_fft
358              indsy=n1*(nd2*k3+ ffti2_local(k2+1) -1)+k1+1        ! this is indsy in the current proc
359              phr=phnons(1,iup,imagn);if (rep==1) phr=phr*symafm_used(jsym) !if rep==2, symafm is already included in phnons
360              phi=phnons(2,iup,imagn);if (rep==1) phi=phi*symafm_used(jsym) !(see irrzg.F90)
361              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)
362              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)
363              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)
364              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)
365              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)
366              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)
367              magxsu1=magxsu1+mxr*phr-mxi*phi;magxsu2=magxsu2+mxi*phr+mxr*phi
368              magysu1=magysu1+myr*phr-myi*phi;magysu2=magysu2+myi*phr+myr*phi
369              magzsu1=magzsu1+mzr*phr-mzi*phi;magzsu2=magzsu2+mzi*phr+mzr*phi
370            end if
371          end do
372          rhosu1_arr(3*izone-2)=magxsu1/dble(nsym_used)
373          rhosu1_arr(3*izone-1)=magysu1/dble(nsym_used)
374          rhosu1_arr(3*izone  )=magzsu1/dble(nsym_used)
375          rhosu2_arr(3*izone-2)=magxsu2/dble(nsym_used)
376          rhosu2_arr(3*izone-1)=magysu2/dble(nsym_used)
377          rhosu2_arr(3*izone  )=magzsu2/dble(nsym_used)
378          numpt=numpt+nup
379        end do
380 
381 !      Reduction in case of FFT parallelization
382        if(mpi_enreg%nproc_fft>1)then
383          spaceComm=mpi_enreg%comm_fft
384          call xmpi_sum(rhosu1_arr,spaceComm,ier)
385          call xmpi_sum(rhosu2_arr,spaceComm,ier)
386        end if
387 
388 !      Now symmetrize the magnetization at equivalent points
389 !      Use Transpose[symrel]
390        numpt=0
391        do izone=1,izone_max
392          rep=irrzon(izone,2,imagn)
393          nup=nsym_used/rep
394          do iup=1,nup
395            ind=irrzon(iup+numpt,1,imagn)
396            j=ind-1;j1=modulo(j,n1);j2=modulo(j/n1,n2);j3=j/(n1*n2)
397            if(fftn2_distrib(j2+1)==me_fft)  then ! this ind is to be treated by me_fft
398              r2=ffti2_local(j2+1) - 1
399              ind=n1*(nd2*j3+r2)+j1+1  ! this is ind in the current proc
400              jsym=isymg(iup+numpt)
401              if (jsym==0) then
402                MSG_ERROR("jsym=0 in symrhg !")
403              end if
404              magxsu1=rhosu1_arr(3*izone-2);magxsu2=rhosu2_arr(3*izone-2)
405              magysu1=rhosu1_arr(3*izone-1);magysu2=rhosu2_arr(3*izone-1)
406              magzsu1=rhosu1_arr(3*izone  );magzsu2=rhosu2_arr(3*izone  )
407              phr=phnons(1,iup,imagn);if (rep==1) phr=phr*symafm_used(jsym) !if rep==2, symafm is already included in phnons
408              phi=phnons(2,iup,imagn);if (rep==1) phi=phi*symafm_used(jsym) !(see irrzg.F90)
409              mxr=symrec_cart(1,1,jsym)*magxsu1+symrec_cart(2,1,jsym)*magysu1+symrec_cart(3,1,jsym)*magzsu1
410              mxi=symrec_cart(1,1,jsym)*magxsu2+symrec_cart(2,1,jsym)*magysu2+symrec_cart(3,1,jsym)*magzsu2
411              myr=symrec_cart(1,2,jsym)*magxsu1+symrec_cart(2,2,jsym)*magysu1+symrec_cart(3,2,jsym)*magzsu1
412              myi=symrec_cart(1,2,jsym)*magxsu2+symrec_cart(2,2,jsym)*magysu2+symrec_cart(3,2,jsym)*magzsu2
413              mzr=symrec_cart(1,3,jsym)*magxsu1+symrec_cart(2,3,jsym)*magysu1+symrec_cart(3,3,jsym)*magzsu1
414              mzi=symrec_cart(1,3,jsym)*magxsu2+symrec_cart(2,3,jsym)*magysu2+symrec_cart(3,3,jsym)*magzsu2
415              magngx(1,ind)=mxr*phr+mxi*phi
416              magngx(2,ind)=mxi*phr-mxr*phi
417              magngy(1,ind)=myr*phr+myi*phi
418              magngy(2,ind)=myi*phr-myr*phi
419              magngz(1,ind)=mzr*phr+mzi*phi
420              magngz(2,ind)=mzi*phr-mzr*phi
421            end if
422          end do
423          numpt=numpt+nup
424        end do
425        ABI_DEALLOCATE(isymg)
426        ABI_DEALLOCATE(rhosu1_arr)
427        ABI_DEALLOCATE(rhosu2_arr)
428        ABI_DEALLOCATE(symrec_cart)
429        ABI_DEALLOCATE(symrel_cart)
430        ABI_DEALLOCATE(symafm_used)
431 
432      end if ! nspden==4
433 
434      call timab(17,2,tsec)
435 
436 !    Pull out full or spin up density, now symmetrized
437      call fourdp(cplex,rhog,work,1,mpi_enreg,nfft,ngfft,paral_kgb,0)
438      rhor(:,ispden)=work(:)
439      if (nspden==4) then
440        call fourdp(cplex,magngx,work,1,mpi_enreg,nfft,ngfft,paral_kgb,0)
441        rhor(:,2)=work(:)
442        call fourdp(cplex,magngy,work,1,mpi_enreg,nfft,ngfft,paral_kgb,0)
443        rhor(:,3)=work(:)
444        call fourdp(cplex,magngz,work,1,mpi_enreg,nfft,ngfft,paral_kgb,0)
445        rhor(:,4)=work(:)
446        ABI_DEALLOCATE(magngx)
447        ABI_DEALLOCATE(magngy)
448        ABI_DEALLOCATE(magngz)
449      end if
450 
451    end do ! ispden
452 
453  end if !  End on the condition nsym==1
454 
455  ABI_DEALLOCATE(work)
456 
457  contains
458 
459    function map_symrhg(j1,n1)
460 
461 
462 !This section has been created automatically by the script Abilint (TD).
463 !Do not modify the following lines by hand.
464 #undef ABI_FUNC
465 #define ABI_FUNC 'map_symrhg'
466 !End of the abilint section
467 
468    integer :: map_symrhg
469    integer,intent(in) :: j1,n1
470    map_symrhg=mod(n1+mod(j1,n1),n1)
471  end function map_symrhg
472 
473 end subroutine symrhg