TABLE OF CONTENTS
ABINIT/symrhg [ 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