TABLE OF CONTENTS


ABINIT/hartredq.F90 [ Functions ]

[ Top ] [ Functions ]

NAME

  hartredq.F90

FUNCTION

  Given rho(G), compute the q-gradient of the Hartree potential at q=0
  (=FFT of -rho(G)*G_qdir/pi**2/|G|**4 ) -> Cartesian coordinates
  The calculation is performed in reduced reciprocal space coordinates.

COPYRIGHT

  Copyright (C) 2021-2022 ABINIT group (FIXME: add author)
  This file is distributed under the terms of the
  GNU General Public License, see ~abinit/COPYING
  or http://www.gnu.org/copyleft/gpl.txt .

INPUTS

  cplex= if 1, vqgradhartr is REAL, if 2, vqgradhartr is COMPLEX
  gmet(3,3)=metrix tensor in G space in Bohr**-2.
  gprimd(3,3)=reciprocal space dimensional primitive translations
  gsqcut=cutoff value on G**2 for sphere inside fft box.
         (gsqcut=(boxcut**2)*ecut/(2.d0*(Pi**2))
  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/input_variables/vargs.htm#ngfft
  qdir= indicates the direction of the q-gradient (1,2 or 3)
  rhog(2,nfft)=electron density in G space

OUTPUT

  vqgradhart(cplex*nfft)=q-gradient of the Hartree potential at q=0in real space, either REAL or COMPLEX

SIDE EFFECTS

NOTES

SOURCE

2647 subroutine hartredq(cplex,gmet,gsqcut,mpi_enreg,nfft,ngfft,qdir,rhog,vqgradhart)
2648 
2649  implicit none
2650 
2651 !Arguments ------------------------------------
2652 !scalars
2653  integer,intent(in) :: cplex,nfft,qdir
2654  real(dp),intent(in) :: gsqcut
2655  type(MPI_type),intent(in) :: mpi_enreg
2656 !arrays
2657  integer,intent(in) :: ngfft(18)
2658  real(dp),intent(in) :: gmet(3,3),rhog(2,nfft)
2659  real(dp),intent(out) :: vqgradhart(cplex*nfft)
2660 
2661 !Local variables-------------------------------
2662 !scalars
2663  integer,parameter :: im=2,re=1
2664  integer :: i1,i2,i23,i2_local,i3
2665  integer :: id1,id2,id3,ig1,ig2,ig3,ii,ii1,me_fft,n1,n2,n3,nproc_fft
2666  real(dp) :: cutoff,gfact,gnorm,num
2667  real(dp), parameter :: piinv2= piinv*two
2668  real(dp),parameter :: tolfix=1.000000001e0_dp
2669 !arrays
2670  integer, ABI_CONTIGUOUS pointer :: fftn2_distrib(:),ffti2_local(:)
2671  integer, ABI_CONTIGUOUS pointer :: fftn3_distrib(:),ffti3_local(:)
2672  real(dp),allocatable :: work1(:,:)
2673  real(dp) :: gvec(3)
2674 
2675 ! *************************************************************************
2676 
2677  DBG_ENTER("COLL")
2678 
2679  n1=ngfft(1); n2=ngfft(2); n3=ngfft(3)
2680  nproc_fft = mpi_enreg%nproc_fft; me_fft = mpi_enreg%me_fft
2681 
2682 !Get the distrib associated with this fft_grid
2683  call ptabs_fourdp(mpi_enreg,n2,n3,fftn2_distrib,ffti2_local,fftn3_distrib,ffti3_local)
2684 
2685 !Initialize a few quantities
2686  cutoff=gsqcut*tolfix
2687  ABI_MALLOC(work1,(2,nfft))
2688  id1=n1/2+2;id2=n2/2+2;id3=n3/2+2
2689 
2690 !Triple loop on each dimension
2691  do i3=1,n3
2692    ig3=i3-(i3/id3)*n3-1
2693 
2694    do i2=1,n2
2695      ig2=i2-(i2/id2)*n2-1
2696 
2697      if (fftn2_distrib(i2) == me_fft) then
2698        i2_local = ffti2_local(i2)
2699        i23=n1*(i2_local-1 +(n2/nproc_fft)*(i3-1))
2700        !Do the test that eliminates the Gamma point outside of the inner loop
2701        ii1=1
2702        if(i23==0 .and. ig2==0 .and. ig3==0)then
2703          ii1=2
2704          work1(re,1+i23)=zero
2705          work1(im,1+i23)=zero
2706        end if
2707 
2708        ! Final inner loop on the first dimension (note the lower limit)
2709        do i1=ii1,n1
2710          ig1=i1-(i1/id1)*n1-1
2711          ii=i1+i23
2712 
2713          gvec=(/ig1,ig2,ig3/)
2714          gnorm=normv(gvec,gmet,'r') !'r' is to avoid the 2pi scalling
2715 
2716          if (gnorm**2<=cutoff) then
2717            num=dot_product(gmet(qdir,:),gvec(:))
2718            gfact=piinv2*num/gnorm**4
2719            work1(re,ii)=-rhog(re,ii)*gfact
2720            work1(im,ii)=-rhog(im,ii)*gfact
2721          else
2722            work1(re,ii)=zero
2723            work1(im,ii)=zero
2724          end if
2725 
2726        end do ! End loop on i1
2727      end if
2728 
2729    end do ! End loop on i2
2730  end do ! End loop on i3
2731 
2732  ! Fourier Transform the q-gradient of the hartree potential, in reciprocal space it was stored in work1
2733  call fourdp(cplex,work1,vqgradhart,1,mpi_enreg,nfft,1,ngfft,0)
2734 
2735  ABI_FREE(work1)
2736 
2737  DBG_EXIT("COLL")
2738 
2739 end subroutine hartredq

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-2022 ABINIT group (XG, BA, MT, DRH, DCA, GMR, MJV, JWZ)
  This file is distributed under the terms of the
  GNU General Public License, see ~abinit/COPYING
  or http://www.gnu.org/copyleft/gpl.txt .

SOURCE

17 #if defined HAVE_CONFIG_H
18 #include "config.h"
19 #endif
20 
21 #include "abi_common.h"
22 
23 module m_spacepar
24 
25  use defs_basis
26  use m_abicore
27  use m_errors
28  use m_xmpi
29  use m_sort
30 
31  use m_time,            only : timab
32  use defs_abitypes,     only : MPI_type
33  use m_symtk,           only : mati3inv, chkgrp, symdet, symatm, matr3inv
34  use m_geometry,        only : metric, normv, symredcart,wedge_basis,wedge_product
35  use m_gtermcutoff,     only : termcutoff
36  use m_mpinfo,          only : ptabs_fourdp
37  use m_fft,             only : zerosym, fourdp
38 
39  implicit none
40 
41  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)

OUTPUT

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

SOURCE

547 subroutine hartre(cplex,gsqcut,icutcoul,izero,mpi_enreg,nfft,ngfft,nkpt,&
548                  &rcut,rhog,rprimd,vcutgeo,vhartr,&
549                  &qpt) ! Optional arguments
550 
551 !Arguments ------------------------------------
552 !scalars
553  integer,intent(in) :: cplex,icutcoul,izero,nfft,nkpt
554  real(dp),intent(in) :: gsqcut,rcut
555  type(MPI_type),intent(in) :: mpi_enreg
556 !arrays
557  integer,intent(in) :: ngfft(18)
558  real(dp),intent(in) :: rprimd(3,3),rhog(2,nfft),vcutgeo(3)
559  real(dp),intent(in),optional :: qpt(3)
560  real(dp),intent(out) :: vhartr(cplex*nfft)
561 
562 !Local variables-------------------------------
563 !scalars
564  integer,parameter :: im=2,re=1
565  integer :: i1,i2,i23,i2_local,i3,id1,id2,id3
566  integer :: ig,ig1min,ig1,ig1max,ig2,ig2min,ig2max,ig3,ig3min,ig3max
567  integer :: ii,ii1,ing,n1,n2,n3,qeq0,qeq05,me_fft,nproc_fft
568  real(dp),parameter :: tolfix=1.000000001e0_dp
569  real(dp) :: cutoff,den,gqg2p3,gqgm12,gqgm13,gqgm23,gs,gs2,gs3,ucvol
570  character(len=500) :: message
571 !arrays
572  integer :: id(3)
573  integer, ABI_CONTIGUOUS pointer :: fftn2_distrib(:),ffti2_local(:)
574  integer, ABI_CONTIGUOUS pointer :: fftn3_distrib(:),ffti3_local(:)
575  real(dp) :: gmet(3,3),gprimd(3,3),qpt_(3),rmet(3,3),tsec(2)
576  real(dp),allocatable :: gcutoff(:)
577  real(dp),allocatable :: gq(:,:),work1(:,:)
578 
579 ! *************************************************************************
580 
581  ! Keep track of total time spent in hartre
582  call timab(10,1,tsec)
583 
584  ! Check that cplex has an allowed value
585  if(cplex/=1 .and. cplex/=2)then
586    write(message, '(a,i0,a,a)' )&
587    'From the calling routine, cplex=',cplex,ch10,&
588    'but the only value allowed are 1 and 2.'
589    ABI_BUG(message)
590  end if
591 
592  call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
593 
594  n1=ngfft(1); n2=ngfft(2); n3=ngfft(3)
595  nproc_fft = mpi_enreg%nproc_fft; me_fft = mpi_enreg%me_fft
596 
597  ! Get the distrib associated with this fft_grid
598  call ptabs_fourdp(mpi_enreg,n2,n3,fftn2_distrib,ffti2_local,fftn3_distrib,ffti3_local)
599 
600  ! Initialize a few quantities
601  cutoff=gsqcut*tolfix
602  if(present(qpt))then
603    qpt_=qpt
604  else
605    qpt_=zero
606  end if
607  qeq0=0
608  if(qpt_(1)**2+qpt_(2)**2+qpt_(3)**2<1.d-15) qeq0=1
609  qeq05=0
610  if (qeq0==0) then
611    if (abs(abs(qpt_(1))-half)<tol12.or.abs(abs(qpt_(2))-half)<tol12.or.abs(abs(qpt_(3))-half)<tol12) qeq05=1
612  end if
613 
614  ! If cplex=1 then qpt_ should be 0 0 0
615  if (cplex==1.and. qeq0/=1) then
616    write(message,'(a,3e12.4,a,a)')&
617    'cplex=1 but qpt=',qpt_,ch10,&
618    'qpt should be 0 0 0.'
619    ABI_BUG(message)
620  end if
621 
622  ! If FFT parallelism then qpt should not be 1/2
623  if (nproc_fft>1.and.qeq05==1) then
624    write(message, '(a,3e12.4,a,a)' )&
625    'FFT parallelism selected but qpt',qpt_,ch10,&
626    'qpt(i) should not be 1/2...'
627    ABI_ERROR(message)
628  end if
629 
630  !Initialize Gcut-off array from m_gtermcutoff
631  !ABI_MALLOC(gcutoff,(ngfft(1)*ngfft(2)*ngfft(3)))
632  call termcutoff(gcutoff,gsqcut,icutcoul,ngfft,nkpt,rcut,rprimd,vcutgeo)
633 
634  ! In order to speed the routine, precompute the components of g+q
635  ! Also check if the booked space was large enough...
636  ABI_MALLOC(gq,(3,max(n1,n2,n3)))
637  do ii=1,3
638    id(ii)=ngfft(ii)/2+2
639    do ing=1,ngfft(ii)
640      ig=ing-(ing/id(ii))*ngfft(ii)-1
641      gq(ii,ing)=ig+qpt_(ii)
642    end do
643  end do
644  ig1max=-1;ig2max=-1;ig3max=-1
645  ig1min=n1;ig2min=n2;ig3min=n3
646 
647  ABI_MALLOC(work1,(2,nfft))
648  id1=n1/2+2;id2=n2/2+2;id3=n3/2+2
649 
650  ! Triple loop on each dimension
651  do i3=1,n3
652    ig3=i3-(i3/id3)*n3-1
653    ! Precompute some products that do not depend on i2 and i1
654    gs3=gq(3,i3)*gq(3,i3)*gmet(3,3)
655    gqgm23=gq(3,i3)*gmet(2,3)*2
656    gqgm13=gq(3,i3)*gmet(1,3)*2
657 
658    do i2=1,n2
659      ig2=i2-(i2/id2)*n2-1
660      if (fftn2_distrib(i2) == me_fft) then
661        gs2=gs3+ gq(2,i2)*(gq(2,i2)*gmet(2,2)+gqgm23)
662        gqgm12=gq(2,i2)*gmet(1,2)*2
663        gqg2p3=gqgm13+gqgm12
664 
665        i2_local = ffti2_local(i2)
666        i23=n1*(i2_local-1 +(n2/nproc_fft)*(i3-1))
667        ! Do the test that eliminates the Gamma point outside of the inner loop
668        ii1=1
669        if(i23==0 .and. qeq0==1  .and. ig2==0 .and. ig3==0)then
670          ii1=2
671          work1(re,1+i23)=zero
672          work1(im,1+i23)=zero
673        end if
674 
675        ! Final inner loop on the first dimension (note the lower limit)
676        do i1=ii1,n1
677          gs=gs2+ gq(1,i1)*(gq(1,i1)*gmet(1,1)+gqg2p3)
678          ii=i1+i23
679 
680          if(gs<=cutoff)then
681            ! Identify min/max indexes (to cancel unbalanced contributions later)
682            ! Count (q+g)-vectors with similar norm
683            if ((qeq05==1).and.(izero==1)) then
684              ig1=i1-(i1/id1)*n1-1
685              ig1max=max(ig1max,ig1); ig1min=min(ig1min,ig1)
686              ig2max=max(ig2max,ig2); ig2min=min(ig2min,ig2)
687              ig3max=max(ig3max,ig3); ig3min=min(ig3min,ig3)
688            end if
689 
690            den=piinv/gs*gcutoff(ii)
691            work1(re,ii)=rhog(re,ii)*den
692            work1(im,ii)=rhog(im,ii)*den
693          else
694            ! gs>cutoff
695            work1(re,ii)=zero
696            work1(im,ii)=zero
697          end if
698 
699        end do ! End loop on i1
700      end if
701    end do ! End loop on i2
702  end do ! End loop on i3
703 
704  ABI_FREE(gq)
705 
706  if (izero==1) then
707    ! Set contribution of unbalanced components to zero
708 
709    if (qeq0==1) then !q=0
710      call zerosym(work1,2,n1,n2,n3,comm_fft=mpi_enreg%comm_fft,distribfft=mpi_enreg%distribfft)
711 
712    else if (qeq05==1) then
713      !q=1/2; this doesn't work in parallel
714      ig1=-1;if (mod(n1,2)==0) ig1=1+n1/2
715      ig2=-1;if (mod(n2,2)==0) ig2=1+n2/2
716      ig3=-1;if (mod(n3,2)==0) ig3=1+n3/2
717      if (abs(abs(qpt_(1))-half)<tol12) then
718        if (abs(ig1min)<abs(ig1max)) ig1=abs(ig1max)
719        if (abs(ig1min)>abs(ig1max)) ig1=n1-abs(ig1min)
720      end if
721      if (abs(abs(qpt_(2))-half)<tol12) then
722        if (abs(ig2min)<abs(ig2max)) ig2=abs(ig2max)
723        if (abs(ig2min)>abs(ig2max)) ig2=n2-abs(ig2min)
724      end if
725      if (abs(abs(qpt_(3))-half)<tol12) then
726        if (abs(ig3min)<abs(ig3max)) ig3=abs(ig3max)
727        if (abs(ig3min)>abs(ig3max)) ig3=n3-abs(ig3min)
728      end if
729      call zerosym(work1,2,n1,n2,n3,ig1=ig1,ig2=ig2,ig3=ig3,&
730        comm_fft=mpi_enreg%comm_fft,distribfft=mpi_enreg%distribfft)
731    end if
732  end if
733 
734  ! Fourier Transform Vhartree. Vh in reciprocal space was stored in work1
735  call fourdp(cplex,work1,vhartr,1,mpi_enreg,nfft,1,ngfft,0)
736 
737  ABI_FREE(gcutoff)
738  ABI_FREE(work1)
739 
740  call timab(10,2,tsec)
741 
742 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,

SOURCE

1298 subroutine hartrestr(gsqcut,idir,ipert,mpi_enreg,natom,nfft,ngfft,rhog,rprimd,vhartr1)
1299 
1300 !Arguments ------------------------------------
1301 !scalars
1302  integer,intent(in) :: idir,ipert,natom,nfft
1303  real(dp),intent(in) :: gsqcut
1304  type(MPI_type),intent(in) :: mpi_enreg
1305 !arrays
1306  integer,intent(in) :: ngfft(18)
1307  real(dp),intent(in) :: rhog(2,nfft),rprimd(3,3)
1308  real(dp),intent(out) :: vhartr1(nfft)
1309 
1310 !Local variables-------------------------------
1311 !scalars
1312  integer,parameter :: im=2,re=1
1313  integer :: i1,i2,i23,i3,id2,id3,ig,ig2,ig3,ii,ii1,ing,istr,ka,kb,n1,n2,n3
1314  real(dp),parameter :: tolfix=1.000000001_dp
1315  real(dp) :: cutoff,ddends,den,dgsds,gqg2p3,gqgm12,gqgm13,gqgm23,gs,gs2,gs3
1316  real(dp) :: term,ucvol
1317  character(len=500) :: message
1318 !arrays
1319  integer,save :: idx(12)=(/1,1,2,2,3,3,3,2,3,1,2,1/)
1320  integer :: id(3)
1321  integer, ABI_CONTIGUOUS pointer :: fftn2_distrib(:),ffti2_local(:)
1322  integer, ABI_CONTIGUOUS pointer :: fftn3_distrib(:),ffti3_local(:)
1323  real(dp) :: dgmetds(3,3),gmet(3,3),gprimd(3,3),gqr(3),rmet(3,3)
1324  real(dp),allocatable :: gq(:,:),work1(:,:)
1325 
1326 ! *************************************************************************
1327 
1328  if( .not. (ipert==natom+3 .or. ipert==natom+4))then
1329    write(message, '(a,i0,a,a)' )&
1330 &   'From the calling routine, ipert=',ipert,ch10,&
1331 &   'so this routine for the strain perturbation should not be called.'
1332    ABI_BUG(message)
1333  end if
1334 
1335  call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
1336 
1337  n1=ngfft(1) ; n2=ngfft(2) ; n3=ngfft(3)
1338 
1339 !Get the distrib associated with this fft_grid
1340  call ptabs_fourdp(mpi_enreg,n2,n3,fftn2_distrib,ffti2_local,fftn3_distrib,ffti3_local)
1341 
1342 !Initialize a few quantities
1343  cutoff=gsqcut*tolfix
1344 
1345  istr=idir + 3*(ipert-natom-3)
1346 
1347  if(istr<1 .or. istr>6)then
1348    write(message, '(a,i10,a,a,a)' )&
1349 &   'Input dir gives istr=',istr,' not allowed.',ch10,&
1350 &   'Possible values are 1,2,3,4,5,6 only.'
1351    ABI_BUG(message)
1352  end if
1353 
1354  ka=idx(2*istr-1);kb=idx(2*istr)
1355  do ii = 1,3
1356    dgmetds(:,ii)=-(gprimd(ka,:)*gprimd(kb,ii)+gprimd(kb,:)*gprimd(ka,ii))
1357  end do
1358 !For historical reasons:
1359  dgmetds(:,:)=0.5_dp*dgmetds(:,:)
1360 
1361 !In order to speed the routine, precompute the components of g+q
1362 !Also check if the booked space was large enough...
1363  ABI_MALLOC(gq,(3,max(n1,n2,n3)))
1364  do ii=1,3
1365    id(ii)=ngfft(ii)/2+2
1366    do ing=1,ngfft(ii)
1367      ig=ing-(ing/id(ii))*ngfft(ii)-1
1368      gq(ii,ing)=ig
1369    end do
1370  end do
1371 
1372  ABI_MALLOC(work1,(2,nfft))
1373  id2=n2/2+2
1374  id3=n3/2+2
1375 !Triple loop on each dimension
1376  do i3=1,n3
1377    ig3=i3-(i3/id3)*n3-1
1378 !  Precompute some products that do not depend on i2 and i1
1379    gqr(3)=gq(3,i3)
1380    gs3=gq(3,i3)*gq(3,i3)*gmet(3,3)
1381    gqgm23=gq(3,i3)*gmet(2,3)*2
1382    gqgm13=gq(3,i3)*gmet(1,3)*2
1383 
1384    do i2=1,n2
1385      if (fftn2_distrib(i2)==mpi_enreg%me_fft) then
1386        gqr(2)=gq(2,i2)
1387        gs2=gs3+ gq(2,i2)*(gq(2,i2)*gmet(2,2)+gqgm23)
1388        gqgm12=gq(2,i2)*gmet(1,2)*2
1389        gqg2p3=gqgm13+gqgm12
1390        ig2=i2-(i2/id2)*n2-1
1391 !      i23=n1*((i2-1)+n2*(i3-1))
1392        i23=n1*((ffti2_local(i2)-1)+(n2/mpi_enreg%nproc_fft)*(i3-1))
1393 !      Do the test that eliminates the Gamma point outside
1394 !      of the inner loop
1395        ii1=1
1396        if(i23==0  .and. ig2==0 .and. ig3==0)then
1397          ii1=2
1398          work1(re,1+i23)=0.0_dp
1399          work1(im,1+i23)=0.0_dp
1400        end if
1401 
1402 !      Final inner loop on the first dimension
1403 !      (note the lower limit)
1404        do i1=ii1,n1
1405          gs=gs2+ gq(1,i1)*(gq(1,i1)*gmet(1,1)+gqg2p3)
1406          ii=i1+i23
1407          if(gs<=cutoff)then
1408            den=piinv/gs
1409            gqr(1)=gq(1,i1)
1410            dgsds=&
1411 &           (gqr(1)*(dgmetds(1,1)*gqr(1)+dgmetds(1,2)*gqr(2)+dgmetds(1,3)*gqr(3))+  &
1412 &           gqr(2)*(dgmetds(2,1)*gqr(1)+dgmetds(2,2)*gqr(2)+dgmetds(2,3)*gqr(3))+  &
1413 &           gqr(3)*(dgmetds(3,1)*gqr(1)+dgmetds(3,2)*gqr(2)+dgmetds(3,3)*gqr(3)) )
1414            ddends=-piinv*dgsds/gs**2
1415            if(istr<=3)then
1416              term=2.0_dp*ddends-den
1417            else
1418              term=2.0_dp*ddends
1419            end if
1420            work1(re,ii)=rhog(re,ii)*term
1421            work1(im,ii)=rhog(im,ii)*term
1422          else
1423            work1(re,ii)=0.0_dp
1424            work1(im,ii)=0.0_dp
1425          end if
1426 
1427        end do ! End loop on i1
1428      end if
1429    end do ! End loop on i2
1430  end do !  End loop on i3
1431 
1432  ABI_FREE(gq)
1433 
1434 !Fourier Transform Vhartree.
1435 !Vh in reciprocal space was stored in work1
1436  call fourdp(1,work1,vhartr1,1,mpi_enreg,nfft,1,ngfft,0)
1437 
1438  ABI_FREE(work1)
1439 
1440 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

SOURCE

2002 subroutine irrzg(irrzon,nspden,nsppol,nsym,n1,n2,n3,phnons,symafm,symrel,tnons)
2003 
2004 !Arguments ------------------------------------
2005 !scalars
2006  integer,intent(in) :: n1,n2,n3,nspden,nsppol,nsym
2007 !arrays
2008  integer,intent(in) :: symafm(nsym),symrel(3,3,nsym)
2009  integer,intent(out) :: irrzon(n1*n2*n3,2,(nspden/nsppol)-3*(nspden/4))
2010  real(dp),intent(in) :: tnons(3,nsym)
2011  real(dp),intent(out) :: phnons(2,n1*n2*n3,(nspden/nsppol)-3*(nspden/4))
2012 
2013 !Local variables-------------------------------
2014 !scalars
2015  integer :: i1,i2,i3,id1,id2,id3,ifft,imagn,ind1,ind2,ipt,irep,isym,izone
2016  integer :: izonemax,j1,j2,j3,jj,k1,k2,k3,l1,l2,l3,nfftot,npt,nsym_used
2017  integer :: nzone,setzer,sppoldbl
2018  real(dp) :: arg,ph1i,ph1r,ph2i,ph2r,tau1,tau2,tau3
2019  logical,parameter :: afm_noncoll=.true. ! TRUE if antiferro symmetries are used in non-collinear magnetism
2020  character(len=500) :: message
2021 !arrays
2022  integer,allocatable :: class(:),iperm(:),symafm_used(:),symrel_used(:,:,:)
2023  integer,allocatable :: work1(:)
2024  real(dp),allocatable :: tnons_used(:,:),work2(:,:)
2025 
2026 ! *************************************************************************
2027 
2028  ABI_MALLOC(class,(nsym))
2029  ABI_MALLOC(iperm,(nsym))
2030  ABI_MALLOC(work1,(n1*n2*n3))
2031  ABI_MALLOC(work2,(2,n1*n2*n3))
2032 
2033  nfftot=n1*n2*n3
2034 
2035  id1=n1/2+2
2036  id2=n2/2+2
2037  id3=n3/2+2
2038 
2039  sppoldbl=nspden/nsppol;if (nspden==4) sppoldbl=1
2040 
2041  do imagn=1,sppoldbl
2042 
2043 !  Treat in a similar way the case of the full group and the non-magnetic subgroup
2044    nsym_used=0
2045    do isym=1,nsym
2046      if( (imagn==1 .and. sppoldbl==2) .or. symafm(isym)==1 .or. &
2047 &     ((nspden==4).and.afm_noncoll) )then
2048        nsym_used=nsym_used+1
2049      end if
2050    end do
2051 
2052    if(imagn==2 .and. nsym_used/=nsym/2)then
2053      write(message, '(a,a,a,a,a,i4,a,i0)' )&
2054 &     '  The number of ferromagnetic symmetry operations must be',ch10,&
2055 &     '  half the total number of operations, while it is observed that',ch10,&
2056 &     '  nsym=',nsym,' and nsym_magn=',nsym_used
2057      ABI_BUG(message)
2058    end if
2059 
2060    ABI_MALLOC(symafm_used,(nsym_used))
2061    ABI_MALLOC(symrel_used,(3,3,nsym_used))
2062    ABI_MALLOC(tnons_used,(3,nsym_used))
2063 
2064    nsym_used=0
2065    do isym=1,nsym
2066      if( (imagn==1 .and. sppoldbl==2) .or. symafm(isym)==1 .or.  &
2067 &     ((nspden==4).and.afm_noncoll) ) then
2068        nsym_used=nsym_used+1
2069        symrel_used(:,:,nsym_used)=symrel(:,:,isym)
2070        tnons_used(:,nsym_used)=tnons(:,isym)
2071        symafm_used(nsym_used)=symafm(isym)
2072      end if
2073    end do
2074    if ((nspden/=4).or.(.not.afm_noncoll)) symafm_used=1
2075 
2076 
2077 !  Zero out work array--later on, a zero entry will mean that
2078 !  a given grid point has not yet been assigned to an ibz point
2079    work1(1:nfftot)=0
2080    irrzon(:,2,imagn)=0
2081 
2082 !  Initialize at G=0 (always in irreducible zone)
2083    nzone=1
2084    irrzon(1,1,imagn)=1
2085    irrzon(1,2,imagn)=nsym_used
2086 !  Set phase exp(2*Pi*I*G dot tnons) for G=0
2087    phnons(1,1,imagn)=one
2088    phnons(2,1,imagn)=zero
2089    npt=1
2090 !  setting work1(1)=1 indicates that first grid point (G=0) is
2091 !  in the iz (irreducible zone)
2092    work1(1)=1
2093 
2094    ind1=0
2095 
2096 !  Loop over reciprocal space grid points:
2097    do i3=1,n3
2098      do i2=1,n2
2099        do i1=1,n1
2100 
2101          ind1=ind1+1
2102 
2103 !        Check to see whether present grid point is equivalent to
2104 !        any previously identified ibz point--if not, a new ibz point
2105 !        has been found
2106 
2107          if (work1(ind1)==0) then
2108 
2109 !          A new class has been found.
2110 
2111 !          Get location of G vector (grid point) centered at 0 0 0
2112            l3=i3-(i3/id3)*n3-1
2113            l2=i2-(i2/id2)*n2-1
2114            l1=i1-(i1/id1)*n1-1
2115 
2116            do isym=1,nsym_used
2117 
2118 !            Get rotated G vector Gj for each symmetry element
2119 !            -- here we use the TRANSPOSE of symrel_used; assuming symrel_used expresses
2120 !            the rotation in real space, the transpose is then appropriate
2121 !            for G space symmetrization (p. 1172d,e of notes, 2 June 1995).
2122              j1=symrel_used(1,1,isym)*l1+&
2123 &             symrel_used(2,1,isym)*l2+symrel_used(3,1,isym)*l3
2124              j2=symrel_used(1,2,isym)*l1+&
2125 &             symrel_used(2,2,isym)*l2+symrel_used(3,2,isym)*l3
2126              j3=symrel_used(1,3,isym)*l1+&
2127 &             symrel_used(2,3,isym)*l2+symrel_used(3,3,isym)*l3
2128 
2129 !            Map into [0,n-1]
2130              k1=mod(n1+mod(j1,n1),n1)
2131              k2=mod(n2+mod(j2,n2),n2)
2132              k3=mod(n3+mod(j3,n3),n3)
2133 
2134 !            Get linear index of rotated point Gj
2135              ind2=1+k1+n1*(k2+n2*k3)
2136 
2137 !            Store info for new class:
2138              class(isym)=ind2
2139              iperm(isym)=isym
2140 
2141 !            Setting work array element to 1 indicates grid point has been
2142 !            identified with iz point
2143              work1(ind2)=1
2144 
2145 !            End of loop on isym
2146            end do
2147 
2148 !          Sort integers into ascending order in each class
2149 !          (this lumps together G vectors with the same linear index, i.e.
2150 !          groups together symmetries which land on the same Gj)
2151            call sort_int(nsym_used,class,iperm)
2152 
2153 !          Check repetition factor (how many identical copies of Gj occur
2154 !          from all symmetries applied to G)
2155            irep=0
2156            do isym=1,nsym_used
2157              if (class(isym)==class(1)) then
2158                irep=irep+1
2159              end if
2160            end do
2161            ipt=nsym_used/irep
2162 
2163 !          Repetition factor must be divisor of nsym_used:
2164            if (nsym_used/=(ipt*irep)) then
2165              write(message, '(a,i5,a,i6,a,a,a,a,a,a)' )&
2166 &             '  irep=',irep,' not a divisor of nsym_used=',nsym_used,ch10,&
2167 &             ' This usually indicates that',&
2168 &             ' the input symmetries do not form a group.',ch10,&
2169 &             ' Action : check the input symmetries carefully do they',&
2170 &             ' form a group ? If they do, there is a code bug.'
2171              ABI_ERROR(message)
2172            end if
2173 
2174 !          Compute phases for any nonsymmorphic symmetries
2175 !          exp(-2*Pi*I*G dot tau(j)) for each symmetry j with
2176 !          (possibly zero) nonsymmorphic translation tau(j)
2177            do jj=1,nsym_used
2178 !            First get nonsymmorphic translation and see if nonzero
2179 !            (iperm grabs the symmetries in the new order after sorting)
2180              isym=iperm(jj)
2181              tau1=tnons_used(1,isym)
2182              tau2=tnons_used(2,isym)
2183              tau3=tnons_used(3,isym)
2184              if (abs(tau1)>tol12.or.abs(tau2)>tol12&
2185 &             .or.abs(tau3)>tol12) then
2186 !              compute exp(-2*Pi*I*G dot tau) using original G
2187                arg=two_pi*(dble(l1)*tau1+dble(l2)*tau2+dble(l3)*tau3)
2188                work2(1,jj)=cos(arg)
2189                work2(2,jj)=-sin(arg)
2190              else
2191                work2(1,jj)=one
2192                work2(2,jj)=zero
2193              end if
2194            end do
2195 
2196 !          All phases arising from symmetries which map to the same
2197 !          G vector must actually be the same because
2198 !          rho(Strans*G)=exp(2*Pi*I*(G) dot tau_S) rho(G)
2199 !          must be satisfied; if exp(2*Pi*I*(G) dot tau_S) can be different
2200 !          for two different symmetries S which both take G to the same St*G,
2201 !          then the related Fourier components rho(St*G) must VANISH.
2202 !          Thus: set "phase" to ZERO here in that case.
2203 !          The G mappings occur in sets of irep members; if irep=1 then
2204 !          all the G are unique.
2205 !          MT 090212:
2206 !          In the case of antiferromagn. symetries, one can have
2207 !          rho(Strans*G)= -exp(2*Pi*I*(G) dot tau_S) rho(G)
2208 !          (look at the minus !)
2209 !          A special treatment is then operated on phons.
2210 !          The later must be consistent with the use of phnons array
2211 !          in symrhg.F90 routine.
2212 !          XG 001108 :
2213 !          Note that there is a tolerance on the
2214 !          accuracy of tnons, especially when they are found from
2215 !          the symmetry finder (with xred that might be a bit inaccurate)
2216            if (irep > 1) then
2217              do jj=1,nsym_used,irep
2218                setzer=0
2219                ph1r=work2(1,jj);ph1i=work2(2,jj)
2220                do j1=jj,jj+irep-1
2221                  ph2r=work2(1,j1);ph2i=work2(2,j1)
2222                  if (((ph2r+ph1r)**2+(ph2i+ph1i)**2) <= tol14) then
2223                    if (setzer/=1) setzer=-1
2224                  else if (((ph2r-ph1r)**2+(ph2i-ph1i)**2) > tol14) then
2225                    setzer=1
2226                  end if
2227                end do
2228 !              Setzer= 0: phnons are all equal
2229 !              Setzer=-1: phnons are equal in absolute value
2230 !              Setzer= 1: some phnons are different
2231                if (setzer/=0) then
2232                  if (setzer==-1) then
2233                    if (afm_noncoll.and.nspden==4) then
2234                      arg=symafm_used(iperm(jj))
2235                      if (all(symafm_used(iperm(jj:jj+irep-1))==arg)) then
2236                        setzer=1
2237                      else
2238                        do j1=jj,jj+irep-1
2239                          work2(:,j1)=work2(:,j1)*dble(symafm_used(iperm(j1)))
2240                        end do
2241                      end if
2242                    else
2243                      setzer=1
2244                    end if
2245                  end if
2246                  if (setzer==1) work2(:,jj:jj+irep-1)=zero
2247                end if
2248              end do
2249 !            Compress data if irep>1:
2250              jj=0
2251              do isym=1,nsym_used,irep
2252                jj=jj+1
2253                class(jj)=class(isym)
2254                work2(1,jj)=work2(1,isym)
2255                work2(2,jj)=work2(2,isym)
2256              end do
2257            end if
2258 
2259 !          Put new unique points into irrzon array:
2260            irrzon(1+npt:ipt+npt,1,imagn)=class(1:ipt)
2261 
2262 !          Put repetition number into irrzon array:
2263            irrzon(1+nzone,2,imagn)=irep
2264 
2265 !          DEBUG
2266 !          write(std_out,'(a,6i7)' )' irrzg : izone,i1,i2,i3,imagn,irrzon(859,2,1)=',&
2267 !          &      1+nzone,i1,i2,i3,imagn,irrzon(859,2,1)
2268 !          ENDDEBUG
2269 
2270 !          Put phases (or 0) in phnons array:
2271            phnons(:,1+npt:ipt+npt,imagn)=work2(:,1:ipt)
2272 
2273 !          Update number of points in irrzon array:
2274 !          (irep must divide evenly into nsym_used !)
2275            npt=npt+ipt
2276 
2277 !          Update number of classes:
2278            nzone=nzone+1
2279 
2280          end if
2281 !
2282 !        End of loop on reciprocal space points, with indices i1, i2, i3
2283        end do
2284      end do
2285    end do
2286 
2287    ABI_SFREE(symafm_used)
2288    ABI_SFREE(symrel_used)
2289    ABI_SFREE(tnons_used)
2290 
2291  end do ! imagn
2292 
2293 !Make sure number of real space points accounted for equals actual number of grid points
2294  if (npt/=n1*n2*n3) then
2295    write(message, '(a,a,a,a,i10,a,i10,a,a,a,a,a,a,a,a,a)' ) ch10,&
2296 &   ' irrzg : ERROR -',ch10,&
2297 &   '  npt=',npt,' and n1*n2*n3=',n1*n2*n3,' are not equal',ch10,&
2298 &   '  This says that the total of all points in the irreducible',&
2299 &   '  sector in real space',ch10,&
2300 &   '  and all symmetrically equivalent',&
2301 &   '  points, npt, does not equal the actual number',ch10,&
2302 &   '  of real space grid points.'
2303    call wrtout(std_out,message,'COLL')
2304    write(message,'(3a)') &
2305 &   ' This may mean that the input symmetries do not form a group',ch10,&
2306 &   ' Action : check input symmetries carefully for errors.'
2307    ABI_ERROR(message)
2308  end if
2309 
2310 !Perform some checks
2311  do imagn=1,sppoldbl
2312 
2313    do ifft=1,nfftot
2314      if (irrzon(ifft,1,imagn)<1.or.irrzon(ifft,1,imagn)>nfftot) then
2315        write(message,'(a,4i0,a,a)')&
2316 &       '  ifft,irrzon(ifft,1,imagn),nfftot,imagn=',ifft,irrzon(ifft,1,imagn),nfftot,imagn,ch10,&
2317 &       '  =>irrzon goes outside acceptable bounds.'
2318        ABI_BUG(message)
2319      end if
2320    end do
2321 
2322    izonemax=0
2323    do izone=1,nfftot
2324 !    Global bounds
2325      if (irrzon(izone,2,imagn)<0.or.irrzon(izone,2,imagn)>(nsym/imagn)) then
2326        write(message, '(a,5i7,a,a)' )&
2327 &       ' izone,nzone,irrzon(izone,2,imagn),nsym,imagn =',izone,nzone,irrzon(izone,2,imagn),nsym,imagn,ch10,&
2328 &       '  =>irrzon goes outside acceptable bounds.'
2329        ABI_BUG(message)
2330      end if
2331 !    Second index only goes up to nzone
2332      if(izonemax==0)then
2333        if (irrzon(izone,2,imagn)==0)izonemax=izone-1
2334      end if
2335      if(izonemax/=0)then
2336        if (irrzon(izone,2,imagn)/=0) then
2337          message = ' beyond izonemax, irrzon(izone,2,imagn) should be zero'
2338          ABI_BUG(message)
2339        end if
2340      end if
2341    end do
2342 
2343  end do ! imagn
2344 
2345  ABI_FREE(class)
2346  ABI_FREE(iperm)
2347  ABI_FREE(work1)
2348  ABI_FREE(work2)
2349 
2350 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
  (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

SOURCE

 982 subroutine laplacian(gprimd,mpi_enreg,nfft,nfunc,ngfft,rdfuncr,&
 983 &  laplacerdfuncr,rdfuncg_out,laplacerdfuncg_out,g2cart_out,rdfuncg_in,g2cart_in)
 984 
 985 !Arguments ------------------------------------
 986 !scalars
 987  integer,intent(in) :: nfft,nfunc
 988  type(MPI_type),intent(in) :: mpi_enreg
 989 !arrays
 990  integer,intent(in) :: ngfft(18)
 991  real(dp),intent(in) :: gprimd(3,3)
 992  real(dp),intent(inout),optional :: laplacerdfuncr(nfft,nfunc)
 993  real(dp),intent(inout),optional,target :: rdfuncr(nfft,nfunc)
 994  real(dp),intent(in),optional,target :: g2cart_in(nfft) !vz_i
 995  real(dp),intent(out),optional,target :: g2cart_out(nfft)  !vz_i
 996  real(dp),intent(out),optional,target :: laplacerdfuncg_out(2,nfft,nfunc)
 997  real(dp),intent(in),optional,target :: rdfuncg_in(2,nfft,nfunc) !vz_i
 998  real(dp),intent(out),optional,target :: rdfuncg_out(2,nfft,nfunc)
 999 
1000 !Local variables-------------------------------
1001 !scalars
1002  integer :: count,i1,i2,i3,id1,id2,id3,ifft,ifunc,ig1,ig2,ig3,ii1,n1,n2
1003  integer :: n3
1004  real(dp) :: b11,b12,b13,b21,b22,b23,b31,b32,b33
1005 !arrays
1006  integer, ABI_CONTIGUOUS pointer :: fftn2_distrib(:),ffti2_local(:)
1007  integer, ABI_CONTIGUOUS pointer :: fftn3_distrib(:),ffti3_local(:)
1008  real(dp),ABI_CONTIGUOUS pointer :: g2cart(:),laplacerdfuncg(:,:,:),rdfuncg(:,:,:)
1009 
1010 ! *************************************************************************
1011 
1012 !Keep local copy of fft dimensions
1013  n1=ngfft(1)
1014  n2=ngfft(2)
1015  n3=ngfft(3)
1016 
1017  if(present(laplacerdfuncg_out)) then
1018    laplacerdfuncg => laplacerdfuncg_out
1019  else
1020    ABI_MALLOC(laplacerdfuncg,(2,nfft,nfunc))
1021  end if
1022 
1023  ! Get the distrib associated with this fft_grid
1024  call ptabs_fourdp(mpi_enreg,n2,n3,fftn2_distrib,ffti2_local,fftn3_distrib,ffti3_local)
1025 
1026 !change the real density rdfuncr on real space on the real density
1027 !rdfuncg in reciprocal space
1028  if(.not.present(rdfuncg_in)) then
1029    if(present(rdfuncg_out)) then
1030      rdfuncg => rdfuncg_out
1031    else
1032      ABI_MALLOC(rdfuncg,(2,nfft,nfunc))
1033    end if
1034    if(present(rdfuncr)) then
1035      do ifunc=1,nfunc
1036        call fourdp(1,rdfuncg(:,:,ifunc),rdfuncr(:,ifunc),-1,mpi_enreg,nfft,1,ngfft,0)
1037      end do
1038    end if
1039  else
1040    rdfuncg => rdfuncg_in
1041  end if
1042 
1043 !apply the laplacian on laplacerdfuncr
1044 !code from /3xc/xcden.F90
1045 !see also src/5common/hatre.F90 and src/5common/moddiel.F90
1046 !Keep local copy of fft dimensions
1047 !Initialize computation of G^2 in cartesian coordinates
1048  if(.not.present(g2cart_in)) then
1049    if(present(g2cart_out)) then
1050      g2cart => g2cart_out
1051    else
1052      ABI_MALLOC(g2cart,(nfft))
1053    end if
1054    id1=int(n1/2)+2
1055    id2=int(n2/2)+2
1056    id3=int(n3/2)+2
1057    count=0
1058    do i3=1,n3
1059      ifft=(i3-1)*n1*(n2/mpi_enreg%nproc_fft)
1060      ig3=i3-int(i3/id3)*n3-1
1061      do i2=1,n2
1062        if (fftn2_distrib(i2)==mpi_enreg%me_fft) then
1063          ig2=i2-int(i2/id2)*n2-1
1064 
1065          ii1=1
1066          do i1=ii1,n1
1067            ig1=i1-int(i1/id1)*n1-1
1068            ifft=ifft+1
1069 
1070            b11=gprimd(1,1)*real(ig1,dp)
1071            b21=gprimd(2,1)*real(ig1,dp)
1072            b31=gprimd(3,1)*real(ig1,dp)
1073            b12=gprimd(1,2)*real(ig2,dp)
1074            b22=gprimd(2,2)*real(ig2,dp)
1075            b32=gprimd(3,2)*real(ig2,dp)
1076            b13=gprimd(1,3)*real(ig3,dp)
1077            b23=gprimd(2,3)*real(ig3,dp)
1078            b33=gprimd(3,3)*real(ig3,dp)
1079 
1080            g2cart(ifft)=( &
1081 &           (b11+b12+b13)**2&
1082 &           +(b21+b22+b23)**2&
1083 &           +(b31+b32+b33)**2&
1084 &           )
1085            do ifunc=1,nfunc
1086 !            compute the laplacian in Fourier space that is * (i x 2pi x G)**2
1087              laplacerdfuncg(1,ifft,ifunc) = -rdfuncg(1,ifft,ifunc)*g2cart(ifft)*two_pi*two_pi
1088              laplacerdfuncg(2,ifft,ifunc) = -rdfuncg(2,ifft,ifunc)*g2cart(ifft)*two_pi*two_pi
1089            end do
1090          end do
1091        end if
1092      end do
1093    end do
1094    if(.not.present(g2cart_out))  then
1095      ABI_FREE(g2cart)
1096    end if
1097  else
1098    g2cart => g2cart_in
1099    do ifunc=1,nfunc
1100      do ifft=1,nfft
1101 !      compute the laplacian in Fourier space that is * (i x 2pi x G)**2
1102        laplacerdfuncg(1,ifft,ifunc) = -rdfuncg(1,ifft,ifunc)*g2cart(ifft)*two_pi*two_pi
1103        laplacerdfuncg(2,ifft,ifunc) = -rdfuncg(2,ifft,ifunc)*g2cart(ifft)*two_pi*two_pi
1104      end do
1105    end do
1106  end if
1107 
1108 !get the result back into real space
1109  if(present(laplacerdfuncr)) then
1110    do ifunc=1,nfunc
1111      call fourdp(1,laplacerdfuncg(:,:,ifunc),laplacerdfuncr(:,ifunc),1,mpi_enreg,nfft,1,ngfft,0)
1112    end do
1113  end if
1114 
1115 !deallocate pointers
1116  if((.not.present(rdfuncg_in)).and.(.not.present(rdfuncg_in)))  then
1117    ABI_FREE(rdfuncg)
1118  end if
1119  if(.not.present(laplacerdfuncg_out))  then
1120    ABI_FREE(laplacerdfuncg)
1121  end if
1122 
1123 end subroutine laplacian

m_spacepar/make_vectornd [ Functions ]

[ Top ] [ m_spacepar ] [ Functions ]

NAME

 make_vectornd

FUNCTION

 For nuclear dipole moments m, compute vector potential A(r) = (m x (r-R))/|r-R|^3
 in r space. This is done by computing A(G) followed by FFT.

NOTES

 This code is copied and modified from m_spacepar/hartre where a very similar loop
 over G is done followed by FFT to real space

INPUTS

OUTPUT

  vectornd(3,nfft)=Vector potential in real space, along Cartesian directions

SOURCE

 83 subroutine make_vectornd(cplex,gsqcut,izero,mpi_enreg,natom,nfft,ngfft,nspden,nucdipmom,&
 84      & rprimd,vectornd,xred)
 85 
 86 !Arguments ------------------------------------
 87 !scalars
 88  integer,intent(in) :: cplex,izero,natom,nfft,nspden
 89  real(dp),intent(in) :: gsqcut
 90  type(MPI_type),intent(in) :: mpi_enreg
 91 !arrays
 92  integer,intent(in) :: ngfft(18)
 93  real(dp),intent(in) :: nucdipmom(3,natom),rprimd(3,3),xred(3,natom)
 94  real(dp),intent(out) :: vectornd(nfft,nspden,3)
 95 
 96 !Local variables-------------------------------
 97  !scalars
 98  integer,parameter :: im=2,re=1
 99  integer :: i1,i2,i2_local,i23,i3,iatom,id1,id2,id3,ig,ig1,ig2,ig3,ig1max,ig2max,ig3max
100  integer :: ig1min,ig2min,ig3min
101  integer :: ii,ii1,ing,me_fft,n1,n2,n3,nd_atom,nd_atom_tot,nproc_fft
102  real(dp),parameter :: tolfix=1.000000001e0_dp
103  real(dp) :: cutoff,gqgm12,gqg2p3,gqgm23,gqgm13,gs2,gs3,gs,phase,ucvol
104  complex(dpc) :: prefac,cgr
105  !arrays
106  integer :: id(3)
107  integer,allocatable :: nd_list(:)
108  integer, ABI_CONTIGUOUS pointer :: fftn2_distrib(:),ffti2_local(:)
109  integer, ABI_CONTIGUOUS pointer :: fftn3_distrib(:),ffti3_local(:)
110  real(dp) :: gmet(3,3),gprimd(3,3),gqred(3),mcgc(3),rmet(3,3)
111  real(dp) :: rgbasis(3,3,3)
112  real(dp),allocatable :: gq(:,:),nd_m(:,:),ndvecr(:),work1(:,:),work2(:,:),work3(:,:)
113 
114 
115 ! *************************************************************************
116 
117  call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
118 
119  ! make list of atoms with nonzero nuclear dipole moments
120  ! in typical applications only 0 or 1 atoms have nonzero dipoles. This
121  ! code shouldn't even be called if all dipoles are zero.
122  nd_atom_tot = 0
123  do iatom = 1, natom
124     if (any(abs(nucdipmom(:,iatom))>tol8)) then
125        nd_atom_tot = nd_atom_tot + 1
126     end if
127  end do
128 
129  ! construct the basis vectors of the generalized cross product
130  ! real space a, b, c (contained in rprimd)
131  ! reciprocal space a*, b*, c* (contained in gprimd)
132  ! for m x G will need a x a*, a x b* etc (9 a^b type basis vectors)
133  call wedge_basis(gprimd,rprimd,rgbasis)
134 
135  ! note that nucdipmom is input as vectors in atomic units referenced
136  ! to cartesian coordinates
137  ABI_MALLOC(nd_list,(nd_atom_tot))
138  ABI_MALLOC(nd_m,(3,nd_atom_tot))
139  nd_atom_tot = 0
140  do iatom = 1, natom
141     if (any(abs(nucdipmom(:,iatom))>tol8)) then
142        nd_atom_tot = nd_atom_tot + 1
143        nd_list(nd_atom_tot) = iatom
144        ! the following expresses the dipole moment components in units of rprimd translations
145        nd_m(:,nd_atom_tot) = MATMUL(TRANSPOSE(gprimd),nucdipmom(:,iatom))
146     end if
147  end do
148 
149  n1=ngfft(1); n2=ngfft(2); n3=ngfft(3)
150  nproc_fft = mpi_enreg%nproc_fft; me_fft = mpi_enreg%me_fft
151 
152  prefac = -four_pi*j_dpc/(ucvol*two_pi)
153 
154  ! Get the distrib associated with this fft_grid
155  call ptabs_fourdp(mpi_enreg,n2,n3,fftn2_distrib,ffti2_local,fftn3_distrib,ffti3_local)
156 
157  ! Initialize a few quantities
158  cutoff=gsqcut*tolfix
159 
160  ! In order to speed the routine, precompute the components of g+q
161  ! Also check if the booked space was large enough...
162  ABI_MALLOC(gq,(3,max(n1,n2,n3)))
163  do ii=1,3
164    id(ii)=ngfft(ii)/2+2
165    do ing=1,ngfft(ii)
166      ig=ing-(ing/id(ii))*ngfft(ii)-1
167      gq(ii,ing)=ig
168    end do
169  end do
170  ig1max=-1;ig2max=-1;ig3max=-1
171  ig1min=n1;ig2min=n2;ig3min=n3
172 
173  ABI_MALLOC(work1,(2,nfft))
174  ABI_MALLOC(work2,(2,nfft))
175  ABI_MALLOC(work3,(2,nfft))
176  work1=zero; work2=zero; work3=zero
177  id1=n1/2+2;id2=n2/2+2;id3=n3/2+2
178 
179  ! Triple loop on each dimension
180  do i3=1,n3
181    ig3=i3-(i3/id3)*n3-1
182    ! Precompute some products that do not depend on i2 and i1
183    gs3=gq(3,i3)*gq(3,i3)*gmet(3,3)
184    gqgm23=gq(3,i3)*gmet(2,3)*2
185    gqgm13=gq(3,i3)*gmet(1,3)*2
186 
187    do i2=1,n2
188      ig2=i2-(i2/id2)*n2-1
189      if (fftn2_distrib(i2) == me_fft) then
190        gs2=gs3+ gq(2,i2)*(gq(2,i2)*gmet(2,2)+gqgm23)
191        gqgm12=gq(2,i2)*gmet(1,2)*2
192        gqg2p3=gqgm13+gqgm12
193 
194        i2_local = ffti2_local(i2)
195        i23=n1*(i2_local-1 +(n2/nproc_fft)*(i3-1))
196        ! Do the test that eliminates the Gamma point outside of the inner loop
197        ii1=1
198        !if(i23==0 .and. ig2==0 .and. ig3==0)then
199        !  ii1=2
200        !  work1(re,1+i23)=zero
201        !  work1(im,1+i23)=zero
202        !end if
203 
204        ! Final inner loop on the first dimension (note the lower limit)
205        do i1=ii1,n1
206           gs=gs2+ gq(1,i1)*(gq(1,i1)*gmet(1,1)+gqg2p3)
207           ig1 = i1 - (i1/id1)*n1 -1
208           ii=i1+i23
209 
210           gqred(1) = gq(1,i1); gqred(2) = gq(2,i2); gqred(3) = gq(3,i3)
211 
212           if( (gs .LE. cutoff) .AND. (gs .gt. tol8) )then
213 
214              do iatom = 1, nd_atom_tot
215                 nd_atom = nd_list(iatom)
216                 phase = -two_pi*DOT_PRODUCT(xred(:,nd_atom),gqred(:))
217                 cgr = cmplx(cos(phase),sin(phase))
218 
219                 ! cross product m x G
220                 call wedge_product(mcgc,nd_m(:,iatom),gqred,rgbasis)
221 
222                 ! express mcgc relative to rprimd translations. This is done because
223                 ! we wish ultimately to apply A.p to |cwavef>; in getghc_nucdip, the
224                 ! p|cwavef> is done in reduced coordinates so do that here too, because
225                 ! r.G has no need of the metric if both terms are in reduced coords
226                 mcgc = MATMUL(TRANSPOSE(gprimd),mcgc)
227 
228                 work1(re,ii) = work1(re,ii) + real(prefac*cgr*mcgc(1)/gs)
229                 work2(re,ii) = work2(re,ii) + real(prefac*cgr*mcgc(2)/gs)
230                 work3(re,ii) = work3(re,ii) + real(prefac*cgr*mcgc(3)/gs)
231 
232                 work1(im,ii) = work1(im,ii) + aimag(prefac*cgr*mcgc(1)/gs)
233                 work2(im,ii) = work2(im,ii) + aimag(prefac*cgr*mcgc(2)/gs)
234                 work3(im,ii) = work3(im,ii) + aimag(prefac*cgr*mcgc(3)/gs)
235 
236              end do
237           else
238              ! gs>cutoff
239              work1(re,ii)=zero
240              work1(im,ii)=zero
241              work2(re,ii)=zero
242              work2(im,ii)=zero
243              work3(re,ii)=zero
244              work3(im,ii)=zero
245           end if
246 
247        end do ! End loop on i1
248      end if
249    end do ! End loop on i2
250  end do ! End loop on i3
251 
252  ABI_FREE(gq)
253  ABI_FREE(nd_list)
254  ABI_FREE(nd_m)
255 
256  if ( izero .EQ. 1 ) then
257    ! Set contribution of unbalanced components to zero
258 
259     call zerosym(work1,2,n1,n2,n3,comm_fft=mpi_enreg%comm_fft,distribfft=mpi_enreg%distribfft)
260     call zerosym(work2,2,n1,n2,n3,comm_fft=mpi_enreg%comm_fft,distribfft=mpi_enreg%distribfft)
261     call zerosym(work3,2,n1,n2,n3,comm_fft=mpi_enreg%comm_fft,distribfft=mpi_enreg%distribfft)
262 
263  end if
264 
265  !note nspden effect--the nuclear vector potential conains no electron spin flip operator,
266  ! so vectornd(:,2,:) = vectornd(:,1,:) and vectornd(:,3:4,:) = zero
267  vectornd = zero
268  ! Fourier Transform
269  ABI_MALLOC(ndvecr,(cplex*nfft))
270  ndvecr=zero
271  call fourdp(cplex,work1,ndvecr,1,mpi_enreg,nfft,1,ngfft,0)
272  vectornd(:,1,1)=ndvecr(:)
273  if (nspden .GE. 2) vectornd(:,2,1) = ndvecr(:)
274  ABI_FREE(work1)
275 
276  ndvecr=zero
277  call fourdp(cplex,work2,ndvecr,1,mpi_enreg,nfft,1,ngfft,0)
278  vectornd(:,1,2) = ndvecr(:)
279  if (nspden .GE. 2) vectornd(:,2,2) = ndvecr(:)
280  ABI_FREE(work2)
281 
282  ndvecr=zero
283  call fourdp(cplex,work3,ndvecr,1,mpi_enreg,nfft,1,ngfft,0)
284  vectornd(:,1,3) = ndvecr(:)
285  if (nspden .GE. 2) vectornd(:,2,3) = ndvecr(:)
286  ABI_FREE(work3)
287  ABI_FREE(ndvecr)
288 
289 end subroutine make_vectornd

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

SOURCE

770 subroutine meanvalue_g(ar,diag,filter,istwf_k,mpi_enreg,npw,nspinor,vect,vect1,use_ndo,ar_im)
771 
772 !Arguments ------------------------------------
773 !scalars
774  integer,intent(in) :: filter,istwf_k,npw,nspinor,use_ndo
775  real(dp),intent(out) :: ar
776  real(dp),intent(out),optional :: ar_im
777  type(MPI_type),intent(in) :: mpi_enreg
778 !arrays
779  real(dp),intent(in) :: diag(npw),vect(2,npw*nspinor)
780  real(dp),intent(in) :: vect1(2,npw*nspinor)
781 
782 !Local variables-------------------------------
783 !scalars
784  integer :: i1,ierr,ipw,jpw,me_g0
785  character(len=500) :: message
786 !arrays
787 
788 ! *************************************************************************
789  me_g0 = mpi_enreg%me_g0
790 
791  DBG_CHECK(ANY(filter==(/0,1/)),"Wrong filter")
792  DBG_CHECK(ANY(nspinor==(/1,2/)),"Wrong nspinor")
793  DBG_CHECK(ANY(istwf_k==(/(ipw,ipw=1,9)/)),"Wrong istwf_k")
794 
795  if(nspinor==2 .and. istwf_k/=1)then
796    write(message,'(a,a,a,i6,a,i6)')&
797 &   'When istwf_k/=1, nspinor must be 1,',ch10,&
798 &   'however, nspinor=',nspinor,', and istwf_k=',istwf_k
799    ABI_BUG(message)
800  end if
801 
802  if(use_ndo==1 .and. (istwf_k==2 .and.me_g0==1)) then
803    ABI_BUG('use_ndo==1, not tested, use istwfk=1')
804  end if
805 
806  ar=zero
807  if(present(ar_im)) ar_im=zero
808 
809 !Normal storage mode
810  if(istwf_k==1)then
811 
812 !  No filter
813    if(filter==0)then
814 !$OMP PARALLEL DO REDUCTION(+:ar)
815      do ipw=1,npw
816        ar=ar+diag(ipw)*(vect(1,ipw)*vect1(1,ipw)+vect(2,ipw)*vect1(2,ipw))
817      end do
818      if(nspinor==2)then
819 !$OMP PARALLEL DO REDUCTION(+:ar) PRIVATE(jpw)
820        do ipw=1+npw,2*npw
821          jpw=ipw-npw
822          ar=ar+diag(jpw)*(vect(1,ipw)*vect1(1,ipw)+vect(2,ipw)*vect1(2,ipw))
823        end do
824      end if
825      if(use_ndo==1)then
826 !$OMP PARALLEL DO REDUCTION(+:ar_im)
827        do ipw=1,npw
828          ar_im=ar_im+diag(ipw)*(vect1(1,ipw)*vect(2,ipw)-vect1(2,ipw)*vect(1,ipw))
829        end do
830        if(nspinor == 2) then
831 !$OMP PARALLEL DO REDUCTION(+:ar_im) PRIVATE(jpw)
832          do ipw=1+npw,2*npw
833            jpw=ipw-npw
834            ar_im=ar_im+diag(jpw)*(vect1(1,ipw)*vect(2,ipw)-vect1(2,ipw)*vect(1,ipw))
835          end do
836        end if
837      end if
838 
839 !    !$OMP PARALLEL DO REDUCTION(+:ar,ar_im)
840 !    do ipw=1,npw
841 !    ar=ar+diag(ipw)*(vect(1,ipw)*vect1(1,ipw)+vect(2,ipw)*vect1(2,ipw))
842 !    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))
843 !    end do
844 !    if(nspinor==2)then
845 !    !$OMP PARALLEL DO PRIVATE(ipw) REDUCTION(+:ar,ar_im)
846 !    do ipw=1+npw,2*npw
847 !    ar=ar+diag(ipw-npw)*(vect(1,ipw)*vect1(1,ipw)+vect(2,ipw)*vect1(2,ipw))
848 !    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))
849 !    end do
850 !    end if
851    else ! will filter
852 
853 !$OMP PARALLEL DO REDUCTION(+:ar)
854      do ipw=1,npw
855        if(diag(ipw)<huge(0.0d0)*1.d-11)then
856          ar=ar+diag(ipw)*(vect(1,ipw)*vect1(1,ipw)+vect(2,ipw)*vect1(2,ipw))
857        end if
858      end do
859      if(nspinor==2)then
860 !$OMP PARALLEL DO REDUCTION(+:ar) PRIVATE(jpw)
861        do ipw=1+npw,2*npw
862          jpw=ipw-npw
863          if(diag(jpw)<huge(0.0d0)*1.d-11)then
864            ar=ar+diag(jpw)*(vect(1,ipw)*vect1(1,ipw)+vect(2,ipw)*vect1(2,ipw))
865          end if
866        end do
867      end if
868      if(use_ndo==1)then
869        if(.not.present(ar_im)) then
870          ABI_BUG("use_ndo true and ar_im not present")
871        end if
872 !$OMP PARALLEL DO REDUCTION(+:ar_im)
873        do ipw=1,npw
874          if(diag(ipw)<huge(0.0d0)*1.d-11)then
875            ar_im=ar_im+diag(ipw)*(vect1(1,ipw)*vect(2,ipw)-vect1(2,ipw)*vect(1,ipw))
876          end if
877        end do
878        if(nspinor == 2) then
879 !$OMP PARALLEL DO REDUCTION(+:ar_im) PRIVATE(jpw)
880          do ipw=1+npw,2*npw
881            jpw=ipw-npw
882            if(diag(jpw)<huge(0.0d0)*1.d-11)then
883              ar_im=ar_im+diag(jpw)*(vect1(1,ipw)*vect(2,ipw)-vect1(2,ipw)*vect(1,ipw))
884            end if
885          end do
886        end if
887      end if
888 
889 
890 !    !$OMP PARALLEL DO PRIVATE(ipw) REDUCTION(+:ar,ar_im)
891 !    do ipw=1,npw
892 !    if(diag(ipw)<huge(0.0d0)*1.d-11)then
893 !    ar=ar+diag(ipw)*(vect(1,ipw)*vect1(1,ipw)+vect(2,ipw)*vect1(2,ipw))
894 !    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))
895 !    end if
896 !    end do
897 !    if(nspinor==2)then
898 !    !$OMP PARALLEL DO PRIVATE(ipw) REDUCTION(+:ar,ar_im)
899 !    do ipw=1+npw,2*npw
900 !    if(diag(ipw-npw)<huge(0.0d0)*1.d-11)then
901 !    ar=ar+diag(ipw-npw)*(vect(1,ipw)*vect1(1,ipw)+vect(2,ipw)*vect1(2,ipw))
902 !    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))
903 !    end if
904 !    end do
905 !    end if ! nspinor==2
906 
907    end if ! filter==0
908 
909  else if(istwf_k>=2)then
910 
911    if(filter==0)then
912      i1=1
913      if(istwf_k==2 .and. me_g0==1)then ! MPIWF need to know which proc has G=0
914        ar=half*diag(1)*vect(1,1)*vect1(1,1) ; i1=2
915      end if
916 
917 !$OMP PARALLEL DO REDUCTION(+:ar)
918      do ipw=i1,npw
919        ar=ar+diag(ipw)*(vect(1,ipw)*vect1(1,ipw)+vect(2,ipw)*vect1(2,ipw))
920      end do
921 
922    else ! filter/=0
923      i1=1
924      if(istwf_k==2 .and. me_g0==1)then
925        if(diag(1)<huge(0.0d0)*1.d-11)then
926          ar=half*diag(1)*vect(1,1)*vect1(1,1) ; i1=2
927        end if
928      end if
929 
930 !$OMP PARALLEL DO REDUCTION(+:ar)
931      do ipw=i1,npw
932        if(diag(ipw)<huge(0.0d0)*1.d-11)then
933          ar=ar+diag(ipw)*(vect(1,ipw)*vect1(1,ipw)+vect(2,ipw)*vect1(2,ipw))
934        end if
935      end do
936    end if ! filter==0
937 
938    ar=two*ar
939 
940  end if ! istwf_k
941 
942 !MPIWF need to make reduction on ar and ai .
943  if(mpi_enreg%paral_kgb==1)then
944    call xmpi_sum(ar,mpi_enreg%comm_bandspinorfft ,ierr)
945    if(present(ar_im))then
946      call xmpi_sum(ar_im,mpi_enreg%comm_bandspinorfft,ierr)
947    end if
948  end if
949 
950 end subroutine meanvalue_g

m_spacepar/mkunitpawspherepot [ Functions ]

[ Top ] [ m_spacepar ] [ Functions ]

NAME

 mkunitpawspherepot

FUNCTION

 Compute "potential" due to a sphere of radius r_paw at one of the ions, of
 strength 1. This is done for testing the completeness of the PAW projectors.

NOTES

INPUTS

OUTPUT

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

SOURCE

309 subroutine mkunitpawspherepot(cplex,gsqcut,izero,mpi_enreg,natom,nfft,ngfft,&
310      & rpaw,rprimd,vunitpawspherepot,xred,&
311                  the_atom) ! Optional arguments
312 
313 !Arguments ------------------------------------
314 !scalars
315  integer,intent(in) :: cplex,izero,natom,nfft
316  integer,intent(in),optional :: the_atom
317  real(dp),intent(in) :: gsqcut,rpaw
318  type(MPI_type),intent(in) :: mpi_enreg
319 !arrays
320  integer,intent(in) :: ngfft(18)
321  real(dp),intent(in) :: rprimd(3,3),xred(3,natom)
322  real(dp),intent(out) :: vunitpawspherepot(cplex*nfft)
323 
324 !Local variables-------------------------------
325 !scalars
326  integer,parameter :: im=2,re=1
327  integer :: i1,i2,i2_local,i23,i3,iatom,id1,id2,id3,ig,ig1,ig2,ig3,ig1max,ig2max,ig3max
328  integer :: ig1min,ig2min,ig3min
329  integer :: ii,ii1,ing,me_fft,n1,n2,n3,nproc_fft,qeq0,qeq05
330  real(dp),parameter :: tolfix=1.000000001e0_dp
331  real(dp) :: cutoff,gqgm12,gqg2p3,gqgm23,gqgm13,gs2,gs3,gs,ogg0,ogr,ogrpre,phgr,phpaw,ucvol
332  character(len=500) :: message
333 !arrays
334  integer :: id(3)
335  integer, ABI_CONTIGUOUS pointer :: fftn2_distrib(:),ffti2_local(:)
336  integer, ABI_CONTIGUOUS pointer :: fftn3_distrib(:),ffti3_local(:)
337  real(dp) :: gmet(3,3),gprimd(3,3),gqred(3),qpt_(3),rmet(3,3)
338  real(dp),allocatable :: gq(:,:),work1(:,:)
339 
340 
341 ! *************************************************************************
342 
343  ! Check that cplex has an allowed value
344  if(cplex/=1 .and. cplex/=2)then
345    write(message, '(a,i0,a,a)' )&
346    'From the calling routine, cplex=',cplex,ch10,&
347    'but the only value allowed are 1 and 2.'
348    ABI_BUG(message)
349  end if
350 
351  call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
352 
353  ogrpre = four_pi/(ucvol*(two_pi**3))
354  ogg0 = four_pi*(rpaw**3)/(three*ucvol)
355 
356  n1=ngfft(1); n2=ngfft(2); n3=ngfft(3)
357  nproc_fft = mpi_enreg%nproc_fft; me_fft = mpi_enreg%me_fft
358 
359  ! Get the distrib associated with this fft_grid
360  call ptabs_fourdp(mpi_enreg,n2,n3,fftn2_distrib,ffti2_local,fftn3_distrib,ffti3_local)
361 
362  ! Initialize a few quantities
363  cutoff=gsqcut*tolfix
364  !carrying code over from hartre with minimal changes, in present case qpt always zero
365  qpt_=zero
366  qeq0=0
367  if(qpt_(1)**2+qpt_(2)**2+qpt_(3)**2<1.d-15) qeq0=1
368  qeq05=0
369  if (qeq0==0) then
370    if (abs(abs(qpt_(1))-half)<tol12.or.abs(abs(qpt_(2))-half)<tol12.or.abs(abs(qpt_(3))-half)<tol12) qeq05=1
371  end if
372 
373  if (present(the_atom)) then
374    iatom = the_atom
375  else
376    iatom = 1
377  end if
378 
379  ! If cplex=1 then qpt_ should be 0 0 0
380  if (cplex==1.and. qeq0/=1) then
381    write(message,'(a,3e12.4,a,a)')&
382    'cplex=1 but qpt=',qpt_,ch10,&
383    'qpt should be 0 0 0.'
384    ABI_BUG(message)
385  end if
386 
387  ! If FFT parallelism then qpt should not be 1/2
388  if (nproc_fft>1.and.qeq05==1) then
389    write(message, '(a,3e12.4,a,a)' )&
390    'FFT parallelism selected but qpt',qpt_,ch10,&
391    'qpt(i) should not be 1/2...'
392    ABI_ERROR(message)
393  end if
394 
395  ! In order to speed the routine, precompute the components of g+q
396  ! Also check if the booked space was large enough...
397  ABI_MALLOC(gq,(3,max(n1,n2,n3)))
398  do ii=1,3
399    id(ii)=ngfft(ii)/2+2
400    do ing=1,ngfft(ii)
401      ig=ing-(ing/id(ii))*ngfft(ii)-1
402      gq(ii,ing)=ig+qpt_(ii)
403    end do
404  end do
405  ig1max=-1;ig2max=-1;ig3max=-1
406  ig1min=n1;ig2min=n2;ig3min=n3
407 
408  ABI_MALLOC(work1,(2,nfft))
409  work1=zero
410  id1=n1/2+2;id2=n2/2+2;id3=n3/2+2
411 
412  ! Triple loop on each dimension
413  do i3=1,n3
414    ig3=i3-(i3/id3)*n3-1
415    ! Precompute some products that do not depend on i2 and i1
416    gs3=gq(3,i3)*gq(3,i3)*gmet(3,3)
417    gqgm23=gq(3,i3)*gmet(2,3)*2
418    gqgm13=gq(3,i3)*gmet(1,3)*2
419 
420    do i2=1,n2
421      ig2=i2-(i2/id2)*n2-1
422      if (fftn2_distrib(i2) == me_fft) then
423        gs2=gs3+ gq(2,i2)*(gq(2,i2)*gmet(2,2)+gqgm23)
424        gqgm12=gq(2,i2)*gmet(1,2)*2
425        gqg2p3=gqgm13+gqgm12
426 
427        i2_local = ffti2_local(i2)
428        i23=n1*(i2_local-1 +(n2/nproc_fft)*(i3-1))
429        ! Do the test that eliminates the Gamma point outside of the inner loop
430        ii1=1
431        if(i23==0 .and. qeq0==1  .and. ig2==0 .and. ig3==0)then
432          ii1=2
433          work1(re,1+i23)=ogg0
434          work1(im,1+i23)=zero
435        end if
436 
437        ! Final inner loop on the first dimension (note the lower limit)
438        do i1=ii1,n1
439          gs=gs2+ gq(1,i1)*(gq(1,i1)*gmet(1,1)+gqg2p3)
440          ii=i1+i23
441 
442          gqred(1) = gq(1,i1); gqred(2) = gq(2,i2); gqred(3) = gq(3,i3)
443 
444          if(gs<=cutoff)then
445 
446 
447            ! Identify min/max indexes (to cancel unbalanced contributions later)
448            ! Count (q+g)-vectors with similar norm
449            if ((qeq05==1).and.(izero==1)) then
450              ig1=i1-(i1/id1)*n1-1
451              ig1max=max(ig1max,ig1); ig1min=min(ig1min,ig1)
452              ig2max=max(ig2max,ig2); ig2min=min(ig2min,ig2)
453              ig3max=max(ig3max,ig3); ig3min=min(ig3min,ig3)
454            end if
455 
456            phpaw=two_pi*rpaw*gs
457            ogr=ogrpre*(sin(phpaw)-phpaw*cos(phpaw))/(gs**3)
458 
459            phgr = -two_pi*DOT_PRODUCT(xred(:,iatom),gqred(:))
460 
461            work1(re,ii)=cos(phgr)*ogr
462            work1(im,ii)=sin(phgr)*ogr
463          else
464            ! gs>cutoff
465            work1(re,ii)=zero
466            work1(im,ii)=zero
467          end if
468 
469        end do ! End loop on i1
470      end if
471    end do ! End loop on i2
472  end do ! End loop on i3
473 
474  ABI_FREE(gq)
475 
476  if (izero==1) then
477    ! Set contribution of unbalanced components to zero
478 
479    if (qeq0==1) then !q=0
480      call zerosym(work1,2,n1,n2,n3,comm_fft=mpi_enreg%comm_fft,distribfft=mpi_enreg%distribfft)
481 
482    else if (qeq05==1) then
483      !q=1/2; this doesn't work in parallel
484      ig1=-1;if (mod(n1,2)==0) ig1=1+n1/2
485      ig2=-1;if (mod(n2,2)==0) ig2=1+n2/2
486      ig3=-1;if (mod(n3,2)==0) ig3=1+n3/2
487      if (abs(abs(qpt_(1))-half)<tol12) then
488        if (abs(ig1min)<abs(ig1max)) ig1=abs(ig1max)
489        if (abs(ig1min)>abs(ig1max)) ig1=n1-abs(ig1min)
490      end if
491      if (abs(abs(qpt_(2))-half)<tol12) then
492        if (abs(ig2min)<abs(ig2max)) ig2=abs(ig2max)
493        if (abs(ig2min)>abs(ig2max)) ig2=n2-abs(ig2min)
494      end if
495      if (abs(abs(qpt_(3))-half)<tol12) then
496        if (abs(ig3min)<abs(ig3max)) ig3=abs(ig3max)
497        if (abs(ig3min)>abs(ig3max)) ig3=n3-abs(ig3min)
498      end if
499      call zerosym(work1,2,n1,n2,n3,ig1=ig1,ig2=ig2,ig3=ig3,&
500        comm_fft=mpi_enreg%comm_fft,distribfft=mpi_enreg%distribfft)
501    end if
502  end if
503 
504  ! Fourier Transform
505  call fourdp(cplex,work1,vunitpawspherepot,1,mpi_enreg,nfft,1,ngfft,0)
506 
507  ABI_FREE(work1)
508 
509 end subroutine mkunitpawspherepot

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)

SOURCE

1154 subroutine redgr(frin,frredgr,mpi_enreg,nfft,ngfft)
1155 
1156 !Arguments ------------------------------------
1157 !scalars
1158  integer,intent(in) :: nfft
1159  type(MPI_type),intent(in) :: mpi_enreg
1160 !arrays
1161  integer,intent(in) :: ngfft(18)
1162  real(dp),intent(in) :: frin(nfft)
1163  real(dp),intent(out) :: frredgr(nfft,3)
1164 
1165 !Local variables-------------------------------
1166 !scalars
1167  integer :: cplex_tmp,i1,i2,i3,id,idir,ifft,ig,ii,ing,n1,n2,n3
1168 !arrays
1169  real(dp),allocatable :: gg(:,:),wkcmpx(:,:),work(:),workgr(:,:)
1170 
1171 ! *************************************************************************
1172 
1173 !Only real arrays are treated
1174  cplex_tmp=1
1175 
1176 !Keep local copy of fft dimensions
1177  n1=ngfft(1) ; n2=ngfft(2) ; n3=ngfft(3)
1178 
1179 !In order to speed the routine, precompute the components of g, including 2pi factor
1180  ABI_MALLOC(gg,(max(n1,n2,n3),3))
1181  do ii=1,3
1182    id=ngfft(ii)/2+2
1183    do ing=1,ngfft(ii)
1184      ig=ing-(ing/id)*ngfft(ii)-1
1185      gg(ing,ii)=two_pi*ig
1186    end do
1187 !  Note that the G <-> -G symmetry must be maintained
1188    if(mod(ngfft(ii),2)==0)gg(ngfft(ii)/2+1,ii)=zero
1189  end do
1190 
1191  ABI_MALLOC(wkcmpx,(2,nfft))
1192  ABI_MALLOC(work,(nfft))
1193  ABI_MALLOC(workgr,(2,nfft))
1194 
1195 !Obtain rho(G) in wkcmpx from input rho(r)
1196  work(:)=frin(:)
1197 
1198  call fourdp(cplex_tmp,wkcmpx,work,-1,mpi_enreg,nfft,1,ngfft,0)
1199 
1200 !Gradient calculation for three reduced components in turn.
1201 !Code duplicated to remove logic from loops.
1202  do idir=1,3
1203    if(idir==1) then
1204 !$OMP PARALLEL DO PRIVATE(ifft)
1205      do i3=1,n3
1206        ifft=(i3-1)*n1*n2
1207        do i2=1,n2
1208          do i1=1,n1
1209            ifft=ifft+1
1210 !          Multiply by i 2pi G(idir)
1211            workgr(2,ifft)= gg(i1,idir)*wkcmpx(1,ifft)
1212            workgr(1,ifft)=-gg(i1,idir)*wkcmpx(2,ifft)
1213          end do
1214        end do
1215      end do
1216    else if(idir==2) then
1217 !$OMP PARALLEL DO PRIVATE(ifft)
1218      do i3=1,n3
1219        ifft=(i3-1)*n1*n2
1220        do i2=1,n2
1221          do i1=1,n1
1222            ifft=ifft+1
1223 !          Multiply by i 2pi G(idir)
1224            workgr(2,ifft)= gg(i2,idir)*wkcmpx(1,ifft)
1225            workgr(1,ifft)=-gg(i2,idir)*wkcmpx(2,ifft)
1226          end do
1227        end do
1228      end do
1229    else
1230 !$OMP PARALLEL DO PRIVATE(ifft)
1231      do i3=1,n3
1232        ifft=(i3-1)*n1*n2
1233        do i2=1,n2
1234          do i1=1,n1
1235            ifft=ifft+1
1236 !          Multiply by i 2pi G(idir)
1237            workgr(2,ifft)= gg(i3,idir)*wkcmpx(1,ifft)
1238            workgr(1,ifft)=-gg(i3,idir)*wkcmpx(2,ifft)
1239          end do
1240        end do
1241      end do
1242    end if !idir
1243 
1244    call fourdp(cplex_tmp,workgr,work,1,mpi_enreg,nfft,1,ngfft,0)
1245 
1246 !$OMP PARALLEL DO
1247    do ifft=1,nfft
1248      frredgr(ifft,idir)=work(ifft)
1249    end do
1250 
1251  end do !idir
1252 
1253  ABI_FREE(gg)
1254  ABI_FREE(wkcmpx)
1255  ABI_FREE(work)
1256  ABI_FREE(workgr)
1257 
1258 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=information 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

SOURCE

2378 subroutine rotate_rho(cplex, itirev, mpi_enreg, nfft, ngfft, nspden, &
2379 &   rhor1, rhog1_eq, rhor1_eq, symrel1, tnon)
2380 
2381 !args
2382  integer,intent(in) :: cplex, nfft, nspden, itirev
2383  integer,intent(in) :: ngfft(18)
2384 
2385  integer, intent(in) :: symrel1(3,3)
2386  real(dp),intent(in) :: tnon(3)
2387  real(dp),intent(inout) :: rhor1(cplex*nfft,nspden)
2388 
2389  real(dp),intent(out) :: rhog1_eq(2,nfft)
2390  real(dp),intent(out) :: rhor1_eq(cplex*nfft,nspden)
2391 
2392  type(MPI_type),intent(in) :: mpi_enreg
2393 
2394 ! local vars
2395  integer :: id1, id2, id3
2396  integer :: n1, n2, n3, nd2
2397  integer :: l1, l2, l3
2398  integer :: i1, i2, i3, ind1, ind2
2399  integer :: j1, j2, j3
2400  integer :: k1, k2, k3
2401  integer :: nproc_fft, ispden, me_fft
2402  real(dp) :: arg
2403  logical :: t_tnon_nonzero
2404 
2405  real(dp) :: phnon1(2)
2406  real(dp), allocatable :: workg(:,:), workg_eq(:,:)
2407  character(len=500) :: message
2408 
2409 ! *************************************************************************
2410 
2411  n1=ngfft(1);n2=ngfft(2);n3=ngfft(3);nproc_fft=ngfft(10);me_fft=ngfft(11);nd2=n2/nproc_fft
2412 
2413  id1=n1/2+2
2414  id2=n2/2+2
2415  id3=n3/2+2
2416 
2417  rhog1_eq = zero
2418  rhor1_eq = zero
2419 
2420  if (itirev == 2) then
2421    write (message,'(3a,9I4,1a)') 'using time reversal. ',ch10,'Symrel1 = ', symrel1, ch10
2422  else
2423    write (message,'(3a,9I4,1a)') 'no time reversal. ',ch10,'Symrel1 = ', symrel1, ch10
2424  end if
2425  !call wrtout(std_out, message, 'COLL')
2426 
2427  t_tnon_nonzero = (any(abs(tnon) > tol12))
2428 
2429 ! eventually, for FFT parallelization
2430 ! call ptabs_fourdp(mpi_enreg,n2,n3,fftn2_distrib,ffti2_local,fftn3_distrib,ffti3_local)
2431 
2432  ABI_MALLOC(workg,(2,nfft))
2433  ABI_MALLOC(workg_eq,(2,nfft))
2434  do ispden = 1, nspden
2435 
2436 ! fft input rhor1 to reciprocal space: uses work* as a buffer
2437    call fourdp(cplex,workg,rhor1(:,ispden),-1,mpi_enreg,nfft,1,ngfft,0)
2438 
2439 ! below taken from irrzg and setsym
2440 !  Loop over reciprocal space grid points:
2441 !  loop over local points in workg, and get back transform from rhog1,
2442 !  which is presumed complete on each proc
2443    ind1=0
2444    do i3=1,n3
2445      do i2=1,n2
2446 !       if(fftn2_distrib(i2)/=me_fft)  cycle ! this ind is not to be treated by me_fft
2447        do i1=1,n1
2448 
2449          ind1=ind1+1
2450 !       r2=ffti2_local(i2+1) - 1
2451 !       ind=n1*(nd2*i3+r2)+i1+1 !this is ind in the current proc
2452 
2453 !      Get location of G vector (grid point) centered at 0 0 0
2454          l1=i1-(i1/id1)*n1-1
2455          l2=i2-(i2/id2)*n2-1
2456          l3=i3-(i3/id3)*n3-1
2457 
2458 !      Get rotated G vector Gj for each symmetry element
2459 !      -- here we use the TRANSPOSE of symrel1; assuming symrel1 expresses
2460 !      the rotation in real space, the transpose is then appropriate
2461 !      for G space symmetrization (p. 1172d,e of notes, 2 June 1995).
2462          j1=symrel1(1,1)*l1+symrel1(2,1)*l2+symrel1(3,1)*l3
2463          j2=symrel1(1,2)*l1+symrel1(2,2)*l2+symrel1(3,2)*l3
2464          j3=symrel1(1,3)*l1+symrel1(2,3)*l2+symrel1(3,3)*l3
2465 
2466 !      Map into [0,n-1] and then add 1 for array index in [1,n]
2467          k1=1+mod(n1+mod(j1,n1),n1)
2468          k2=1+mod(n2+mod(j2,n2),n2)
2469          k3=1+mod(n3+mod(j3,n3),n3)
2470 
2471 !      Get linear index of rotated point Gj
2472          ind2=k1+n1*((k2-1)+n2*(k3-1))
2473 !       r2=ffti2_local(j2+1) - 1
2474 !       ind=n1*(nd2*j3+r2)+j1+1 !this is ind may be in another proc!!
2475 
2476          phnon1(1) = one
2477          phnon1(2) = zero
2478          if (t_tnon_nonzero) then
2479 !        compute exp(-2*Pi*I*G dot tau) using original G
2480 ! NB: this phase is same as that in irrzg and phnons1, and corresponds to complex conjugate of phase from G to Gj;
2481 ! we use it immediately below, to go _to_ workg(ind1)
2482 ! TODO : replace this with complex powers of exp(2pi tnon(1)) etc...
2483            arg=two_pi*(dble(l1)*tnon(1)+dble(l2)*tnon(2)+dble(l3)*tnon(3))
2484            phnon1(1) = cos(arg)
2485            phnon1(2) =-sin(arg)
2486          end if
2487 
2488 !      rho(Strans*G)=exp(2*Pi*I*(G) dot tau_S) rho(G)
2489          workg_eq (1, ind1) = phnon1(1) * workg(1, ind2) &
2490 &         - phnon1(2) * workg(2, ind2)
2491          workg_eq (2, ind1) = phnon1(1) * workg(2, ind2) &
2492 &         + phnon1(2) * workg(1, ind2)
2493 
2494        end do
2495      end do
2496    end do
2497 
2498 ! accumulate rhog1_eq
2499    if (ispden == 1) rhog1_eq = workg_eq
2500 
2501 ! FFT back to real space to get rhor1_eq
2502 !    Pull out full or spin up density, now symmetrized
2503    call fourdp(cplex,workg_eq,rhor1_eq(:,ispden),1,mpi_enreg,nfft,1,ngfft,0)
2504 
2505  end do !nspden
2506 
2507  ABI_FREE(workg)
2508  ABI_FREE(workg_eq)
2509 
2510 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

SOURCE

2552 subroutine setsym(indsym,irrzon,iscf,natom,nfft,ngfft,nspden,nsppol,nsym,phnons,&
2553 & symafm,symrec,symrel,tnons,typat,xred)
2554 
2555 !Arguments ------------------------------------
2556 !scalars
2557  integer,intent(in) :: iscf,natom,nfft,nspden,nsppol,nsym
2558 !arrays
2559  integer,intent(in) :: ngfft(18),symafm(nsym),symrel(3,3,nsym),typat(natom)
2560  integer,intent(out) :: indsym(4,nsym,natom)
2561  integer,intent(inout) :: irrzon(nfft,2,(nspden/nsppol)-3*(nspden/4)) !vz_i
2562  integer,intent(out) :: symrec(3,3,nsym)
2563  real(dp),intent(in) :: tnons(3,nsym),xred(3,natom)
2564  real(dp),intent(out) :: phnons(2,nfft,(nspden/nsppol)-3*(nspden/4))
2565 
2566 !Local variables-------------------------------
2567 !scalars
2568  integer :: isym,ierr
2569  real(dp) :: tolsym8
2570 !arrays
2571  integer,allocatable :: determinant(:)
2572  real(dp) :: tsec(2)
2573 
2574 ! *************************************************************************
2575 
2576  call timab(6,1,tsec)
2577 
2578 !Check that symmetries have unity determinant
2579  ABI_MALLOC(determinant,(nsym))
2580  call symdet(determinant,nsym,symrel)
2581  ABI_FREE(determinant)
2582 
2583 
2584 !Get the symmetry matrices in terms of reciprocal basis
2585  do isym=1,nsym
2586    call mati3inv(symrel(:,:,isym),symrec(:,:,isym))
2587  end do
2588 
2589 !Check for group closure
2590  call chkgrp(nsym,symafm,symrel,ierr)
2591  ABI_CHECK(ierr==0,"Error in group closure")
2592 
2593  call chkgrp(nsym,symafm,symrec,ierr)
2594  ABI_CHECK(ierr==0,"Error in group closure")
2595 
2596 !Obtain a list of rotated atom labels:
2597  tolsym8=tol8
2598  call symatm(indsym,natom,nsym,symrec,tnons,tolsym8,typat,xred,print_indsym=10)
2599 
2600 !If non-SCF calculation, or nsym==1, do not need IBZ data
2601  if ( (iscf>0 .or. iscf==-3) .and. nsym>1 ) then
2602 !  Locate irreducible zone in reciprocal space for symmetrization:
2603    call irrzg(irrzon,nspden,nsppol,nsym,ngfft(1),ngfft(2),ngfft(3),phnons,symafm,symrel,tnons)
2604  end if
2605 
2606  call timab(6,2,tsec)
2607 
2608 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)
 tnons(3,nsym)=reduced nonsymmorphic translations

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

SOURCE

1487 subroutine symrhg(cplex,gprimd,irrzon,mpi_enreg,nfft,nfftot,ngfft,nspden,nsppol,nsym,&
1488 &                 phnons,rhog,rhor,rprimd,symafm,symrel,tnons)
1489 
1490 !Arguments ------------------------------------
1491 !scalars
1492  integer,intent(in) :: cplex,nfft,nfftot,nspden,nsppol,nsym
1493  type(MPI_type),intent(in) :: mpi_enreg
1494 !arrays
1495  integer,intent(in) :: irrzon(nfftot**(1-1/nsym),2,(nspden/nsppol)-3*(nspden/4)),ngfft(18)
1496  integer,intent(in) :: symafm(nsym),symrel(3,3,nsym)
1497  real(dp),intent(in) :: gprimd(3,3),phnons(2,nfftot**(1-1/nsym),(nspden/nsppol)-3*(nspden/4)),rprimd(3,3)
1498  real(dp),intent(inout) :: rhor(cplex*nfft,nspden)
1499  real(dp),intent(out) :: rhog(2,nfft)
1500  real(dp),intent(in) :: tnons(3,nsym)
1501 
1502 !Local variables-------------------------------
1503 !scalars
1504  integer :: id1,id2,id3,ier,imagn,ind,ind2,indsy,ispden,isym,iup,izone,izone_max,j,j1,j2,j3,jsym
1505  integer :: k1,k2,k3,l1,l2,l3,me_fft
1506  integer :: n1,n2,n3,nd2,nproc_fft,nspden_eff,nsym_used,numpt,nup
1507  integer :: r2,rep,spaceComm
1508  logical,parameter :: afm_noncoll=.true.  ! TRUE if antiferro symmetries are used in non-collinear magnetism
1509  real(dp) :: arg,tau1,tau2,tau3
1510  real(dp) :: magxsu1,magxsu2,magysu1,magysu2,magzsu1,magzsu2,mxi,mxr,myi,myr,mzi,mzr,phi,phr,rhosu1,rhosu2
1511  !character(len=500) :: message
1512 !arrays
1513  integer,allocatable :: isymg(:)
1514  integer, ABI_CONTIGUOUS pointer :: fftn2_distrib(:),ffti2_local(:)
1515  integer, ABI_CONTIGUOUS pointer :: fftn3_distrib(:),ffti3_local(:)
1516  real(dp) :: tsec(2)
1517  real(dp),allocatable :: magngx(:,:),magngy(:,:),magngz(:,:)
1518  real(dp),allocatable :: rhosu1_arr(:),rhosu2_arr(:),work(:)
1519  real(dp),allocatable :: symafm_used(:),symrec_cart(:,:,:),symrel_cart(:,:,:),tnons_used(:,:)
1520 
1521 !*************************************************************************
1522 !
1523 !Note the timing channel 17 excludes the different Fourier transforms
1524 
1525  ABI_MALLOC(work,(cplex*nfft))
1526 
1527 !Special treatment for spin-polarized case
1528  if(nspden==2 .and. nsppol==2) then
1529 !  When nspden=2 and nsppol=2, put total density in first half
1530 !  of rhor array and spin up in second half  (up,down) => (up+down,up)
1531    call timab(17,1,tsec)
1532    work(:)=rhor(:,1)               ! up => work
1533    rhor(:,1)=rhor(:,1)+rhor(:,2)   ! up+down
1534    rhor(:,2)=work(:)               ! work => up
1535    call timab(17,2,tsec)
1536  end if
1537 
1538 !Special treatment for antiferromagnetism case
1539  if(nspden==2 .and. nsppol==1) then
1540    call timab(17,1,tsec)
1541 !  When nspden=2 and nsppol=1, (2*up) => (2*up,up)
1542 !  Indeed, what was delivered to the present routine is a "total" density,
1543 !  obtained from occupation numbers varying between 0 and 2,
1544 !  but for spin up only potential.
1545    rhor(:,2)=half*rhor(:,1)
1546    call timab(17,2,tsec)
1547  end if
1548 
1549 !Special treatment for non-collinear magnetism case
1550  if(nspden==4) then
1551    call timab(17,1,tsec)
1552 !FR the half factors missing are recovered in dfpt_mkvxc_noncoll and dfpt_accrho
1553    rhor(:,1)=rhor(:,1)+rhor(:,4)     !nup+ndown
1554    rhor(:,2)=rhor(:,2)-rhor(:,1)     !mx (n+mx-n)
1555    rhor(:,3)=rhor(:,3)-rhor(:,1)     !my (n+my-n)
1556    rhor(:,4)=rhor(:,1)-two*rhor(:,4) !mz=n-2ndown
1557    call timab(17,2,tsec)
1558  end if
1559 
1560 
1561  if(nsym==1)then
1562 
1563    if(nspden==2 .and. nsppol==1) then ! There must be at least one anti-ferromagnetic operation
1564      ABI_BUG('In the antiferromagnetic case, nsym cannot be 1')
1565    end if
1566 
1567 !  If not using symmetry, still want total density in G space rho(G).
1568 !  Fourier transform (incl normalization) to get rho(G)
1569    work(:)=rhor(:,1)
1570    call fourdp(cplex,rhog,work,-1,mpi_enreg,nfft,1,ngfft,0)
1571  else
1572 
1573 !  Treat either full density, spin-up density or magnetization
1574 !  Note the decrease of ispden to the value 1, in order to finish
1575 !  with rhog of the total density (and not the spin-up density or magnetization)
1576    nspden_eff=nspden;if (nspden==4) nspden_eff=1
1577    do ispden=nspden_eff,1,-1
1578 
1579 !    Prepare the density to be symmetrized, in the reciprocal space
1580      if(nspden==1 .or. nsppol==2 .or. (nspden==4.and.(.not.afm_noncoll)))then
1581        imagn=1
1582        nsym_used=0
1583        do isym=1,nsym
1584          if(symafm(isym)==1)nsym_used=nsym_used+1
1585 !        DEBUG
1586 !        write(std_out,*)' symrhg : isym,symafm(isym)',isym,symafm(isym)
1587 !        ENDDEBUG
1588        end do
1589      else if(nspden==2 .and. nsppol==1)then   ! antiferromagnetic case
1590        imagn=ispden
1591        nsym_used=nsym/ispden
1592      else if (nspden==4) then
1593        imagn=1
1594        nsym_used=nsym/ispden
1595      end if
1596 
1597 !    write(std_out,*)' symrhg : nsym_used=',nsym_used
1598 
1599 !    rhor -fft-> rhog    (rhog is used as work space)
1600 !    Note : it should be possible to reuse rhog in the antiferromagnetic case this would avoid one FFT
1601      work(:)=rhor(:,ispden)
1602      call fourdp(cplex,rhog,work,-1,mpi_enreg,nfft,1,ngfft,0)
1603      if (nspden==4) then
1604        ABI_MALLOC(magngx,(2,nfft))
1605        ABI_MALLOC(magngy,(2,nfft))
1606        ABI_MALLOC(magngz,(2,nfft))
1607        work(:)=rhor(:,2)
1608        call fourdp(cplex,magngx,work,-1,mpi_enreg,nfft,1,ngfft,0)
1609        work(:)=rhor(:,3)
1610        call fourdp(cplex,magngy,work,-1,mpi_enreg,nfft,1,ngfft,0)
1611        work(:)=rhor(:,4)
1612        call fourdp(cplex,magngz,work,-1,mpi_enreg,nfft,1,ngfft,0)
1613      end if
1614 
1615 !    Begins the timing here only , to exclude FFTs
1616      call timab(17,1,tsec)
1617 
1618      n1=ngfft(1);n2=ngfft(2);n3=ngfft(3);nproc_fft=ngfft(10);me_fft=ngfft(11);nd2=n2/nproc_fft
1619 
1620 !    Get the distrib associated with this fft_grid
1621      call ptabs_fourdp(mpi_enreg,n2,n3,fftn2_distrib,ffti2_local,fftn3_distrib,ffti3_local)
1622 
1623 !    The following is only valid for total, up or dn density
1624 !    -------------------------------------------------------
1625 
1626 !    Get maxvalue of izone
1627      izone_max=count(irrzon(:,2,imagn)>0)
1628      ABI_MALLOC(rhosu1_arr,(izone_max))
1629      ABI_MALLOC(rhosu2_arr,(izone_max))
1630 
1631      numpt=0
1632      do izone=1,nfftot
1633 
1634 !      Get repetition number
1635        rep=irrzon(izone,2,imagn)
1636        if(rep==0)exit
1637 
1638 !      Compute number of unique points in this symm class:
1639        nup=nsym_used/rep
1640 
1641 !      Accumulate charge over equivalent points
1642        rhosu1=zero
1643        rhosu2=zero
1644        do iup=1,nup
1645          ind=irrzon(iup+numpt,1,imagn)
1646          j=ind-1;j1=modulo(j,n1);j2=modulo(j/n1,n2);j3=j/(n1*n2);
1647          if(fftn2_distrib(j2+1)==me_fft)  then ! this ind is to be treated by me_fft
1648            r2=ffti2_local(j2+1) - 1
1649            ind=n1*(nd2*j3+r2)+j1+1 !this is ind in the current proc
1650            rhosu1=rhosu1+rhog(1,ind)*phnons(1,iup+numpt,imagn)&
1651 &           -rhog(2,ind)*phnons(2,iup+numpt,imagn)
1652            rhosu2=rhosu2+rhog(2,ind)*phnons(1,iup+numpt,imagn)&
1653 &           +rhog(1,ind)*phnons(2,iup+numpt,imagn)
1654          end if
1655 
1656        end do
1657        rhosu1=rhosu1/dble(nup)
1658        rhosu2=rhosu2/dble(nup)
1659        rhosu1_arr(izone)=rhosu1
1660        rhosu2_arr(izone)=rhosu2
1661 !      Keep index of how many points have been considered:
1662        numpt=numpt+nup
1663 
1664      end do  ! End loop over izone
1665 
1666 !    Reduction in case of FFT parallelization
1667      if(mpi_enreg%nproc_fft>1)then
1668        spaceComm=mpi_enreg%comm_fft
1669        call xmpi_sum(rhosu1_arr,spaceComm,ier)
1670        call xmpi_sum(rhosu2_arr,spaceComm,ier)
1671      end if
1672 
1673 !    Now symmetrize the density
1674      numpt=0
1675      do izone=1,nfftot
1676 
1677 !      Get repetition number
1678        rep=irrzon(izone,2,imagn)
1679        if(rep==0)exit
1680 
1681 !      Compute number of unique points in this symm class:
1682        nup=nsym_used/rep
1683 
1684 !      Define symmetrized rho(G) at equivalent points:
1685        do iup=1,nup
1686          ind=irrzon(iup+numpt,1,imagn)
1687 !        decompose ind-1=n1(n2 j3+ j2)+j1
1688          j=ind-1;j1=modulo(j,n1);j2=modulo(j/n1,n2);j3=j/(n1*n2);
1689          if(fftn2_distrib(j2+1)==me_fft)  then ! this ind is to be treated by me_fft
1690            r2=ffti2_local(j2+1) - 1
1691 !          ind in the proc ind-1=n1(nd2 j3+ r2)+j1
1692            ind=n1*(nd2*j3+r2)+j1+1 !this is ind in the current proc
1693            rhog(1,ind)=rhosu1_arr(izone)*phnons(1,iup+numpt,imagn)&
1694 &           +rhosu2_arr(izone)*phnons(2,iup+numpt,imagn)
1695            rhog(2,ind)=rhosu2_arr(izone)*phnons(1,iup+numpt,imagn)&
1696 &           -rhosu1_arr(izone)*phnons(2,iup+numpt,imagn)
1697          end if
1698        end do
1699 
1700 !      Keep index of how many points have been considered:
1701        numpt=numpt+nup
1702 
1703      end do ! End loop over izone
1704 
1705      ABI_FREE(rhosu1_arr)
1706      ABI_FREE(rhosu2_arr)
1707 
1708 !    The following is only valid for magnetization
1709 !    ---------------------------------------------
1710      if (nspden==4) then
1711 
1712        id1=n1/2+2
1713        id2=n2/2+2
1714        id3=n3/2+2
1715 
1716 !      Transfer symmetries in cartesian coordinates
1717 !      Compute symmetries in reciprocal space in cartesian coordinates
1718        ABI_MALLOC(symrec_cart,(3,3,nsym_used))
1719        ABI_MALLOC(symrel_cart,(3,3,nsym_used))
1720        ABI_MALLOC(symafm_used,(nsym_used))
1721        ABI_MALLOC(tnons_used,(3,nsym_used))
1722        jsym=0
1723        do isym=1,nsym
1724          if (symafm(isym)/=1.and.(.not.afm_noncoll)) cycle
1725          jsym=jsym+1
1726          tnons_used(:,jsym)=tnons(:,isym)
1727          symafm_used(jsym)=dble(symafm(isym))
1728          call symredcart(rprimd,gprimd,symrel_cart(:,:,jsym),symrel(:,:,isym))
1729          call matr3inv(symrel_cart(:,:,jsym),symrec_cart(:,:,jsym))
1730        end do
1731 
1732        numpt=count(irrzon(:,1,imagn)>0)
1733        ABI_MALLOC(isymg,(numpt))
1734        isymg=0
1735        ABI_MALLOC(rhosu1_arr,(3*izone_max))
1736        ABI_MALLOC(rhosu2_arr,(3*izone_max))
1737 
1738 !      Accumulate magnetization over equivalent points
1739 !      Use all symmetries (not only those linking different g points)
1740 !      Use Inverse[Transpose[symrel]]=symrec
1741        numpt=0
1742        do izone=1,izone_max
1743          magxsu1=zero;magxsu2=zero
1744          magysu1=zero;magysu2=zero
1745          magzsu1=zero;magzsu2=zero
1746          ind=irrzon(1+numpt,1,1)
1747          rep=irrzon(izone,2,1)
1748          nup=nsym_used/rep
1749 !        Get coordinates in the range [0,n-1]
1750          j=ind-1;l1=modulo(j,n1);l2=modulo(j/n1,n2);l3=j/(n1*n2)
1751 !        Get location of G vector (grid point) centered at 0 0 0
1752 !TO BE UNCOMMENTED
1753          l3=l3-(l3/id3)*n3
1754          l2=l2-(l2/id2)*n2
1755          l1=l1-(l1/id1)*n1
1756 
1757          jsym=0
1758          do isym=1,nsym
1759            if (symafm(isym)/=1.and.(.not.afm_noncoll)) cycle
1760            jsym=jsym+1
1761 !          The G vectors should transform as vectors in reciprocal space
1762 !          However, one acts with the INVERSE of the symmetry operation => Inverse[symrec]=Transpose[symrel]
1763            j1=symrel(1,1,isym)*l1+symrel(2,1,isym)*l2+symrel(3,1,isym)*l3
1764            j2=symrel(1,2,isym)*l1+symrel(2,2,isym)*l2+symrel(3,2,isym)*l3
1765            j3=symrel(1,3,isym)*l1+symrel(2,3,isym)*l2+symrel(3,3,isym)*l3
1766            k1=map_symrhg(j1,n1);k2=map_symrhg(j2,n2);k3=map_symrhg(j3,n3)
1767            indsy=1+k1+n1*(k2+n2*k3)
1768            ind2=-1;iup=numpt
1769            do while (ind2/=indsy.and.iup<numpt+nup)
1770              iup=iup+1;ind2=irrzon(iup,1,1)
1771            end do
1772            if (ind2/=indsy) then
1773              ABI_ERROR("ind2/=indsy in symrhg !")
1774            end if
1775            if (isymg(iup)==0) isymg(iup)=jsym
1776            if(fftn2_distrib(modulo((indsy-1)/n1,n2) + 1) == me_fft ) then  ! this is indsy is to be treated by me_fft
1777              indsy=n1*(nd2*k3+ ffti2_local(k2+1) -1)+k1+1        ! this is indsy in the current proc
1778 
1779 !            Working on this: the present coding will be detrimental for speed ! cos and sin are recomputed many times !
1780              tau1=tnons_used(1,jsym)
1781              tau2=tnons_used(2,jsym)
1782              tau3=tnons_used(3,jsym)
1783              if (abs(tau1)>tol12.or.abs(tau2)>tol12.or.abs(tau3)>tol12) then
1784 !              Compute exp(-2*Pi*I*G dot tau) using original G (equivalent of phnons in the collinear case)
1785                arg=two_pi*(dble(l1)*tau1+dble(l2)*tau2+dble(l3)*tau3)
1786                phr=cos(arg)
1787                phi=-sin(arg)
1788              else
1789                phr=one
1790                phi=zero
1791              end if
1792              phr=phr*symafm_used(jsym)
1793              phi=phi*symafm_used(jsym)
1794 !TO BE COMMENTED
1795 !            phr=phnons(1,iup,imagn);if (rep==1) phr=phr*symafm_used(jsym) !if rep==2, symafm is already included in phnons
1796 !            phi=phnons(2,iup,imagn);if (rep==1) phi=phi*symafm_used(jsym) !(see irrzg.F90)
1797 
1798 !            The magnetization should transform as a vector in real space
1799 !            However, one acts with the INVERSE of the symmetry operation.
1800 !            => Inverse[symrel_cart] = Transpose[symrel_cart] because symrel_cart is unitary   ?!?!?
1801              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)
1802              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)
1803              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)
1804              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)
1805              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)
1806              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)
1807 
1808 !            mxr=symrel_cart(1,1,jsym)*magngx(1,indsy)+symrel_cart(2,1,jsym)*magngy(1,indsy)+symrel_cart(3,1,jsym)*magngz(1,indsy)
1809 !            mxi=symrel_cart(1,1,jsym)*magngx(2,indsy)+symrel_cart(2,1,jsym)*magngy(2,indsy)+symrel_cart(3,1,jsym)*magngz(2,indsy)
1810 !            myr=symrel_cart(1,2,jsym)*magngx(1,indsy)+symrel_cart(2,2,jsym)*magngy(1,indsy)+symrel_cart(3,2,jsym)*magngz(1,indsy)
1811 !            myi=symrel_cart(1,2,jsym)*magngx(2,indsy)+symrel_cart(2,2,jsym)*magngy(2,indsy)+symrel_cart(3,2,jsym)*magngz(2,indsy)
1812 !            mzr=symrel_cart(1,3,jsym)*magngx(1,indsy)+symrel_cart(2,3,jsym)*magngy(1,indsy)+symrel_cart(3,3,jsym)*magngz(1,indsy)
1813 !            mzi=symrel_cart(1,3,jsym)*magngx(2,indsy)+symrel_cart(2,3,jsym)*magngy(2,indsy)+symrel_cart(3,3,jsym)*magngz(2,indsy)
1814 
1815              magxsu1=magxsu1+mxr*phr-mxi*phi;magxsu2=magxsu2+mxi*phr+mxr*phi
1816              magysu1=magysu1+myr*phr-myi*phi;magysu2=magysu2+myi*phr+myr*phi
1817              magzsu1=magzsu1+mzr*phr-mzi*phi;magzsu2=magzsu2+mzi*phr+mzr*phi
1818            end if
1819          end do
1820          rhosu1_arr(3*izone-2)=magxsu1/dble(nsym_used)
1821          rhosu1_arr(3*izone-1)=magysu1/dble(nsym_used)
1822          rhosu1_arr(3*izone  )=magzsu1/dble(nsym_used)
1823          rhosu2_arr(3*izone-2)=magxsu2/dble(nsym_used)
1824          rhosu2_arr(3*izone-1)=magysu2/dble(nsym_used)
1825          rhosu2_arr(3*izone  )=magzsu2/dble(nsym_used)
1826          numpt=numpt+nup
1827        end do
1828 
1829 !      Reduction in case of FFT parallelization
1830        if(mpi_enreg%nproc_fft>1)then
1831          spaceComm=mpi_enreg%comm_fft
1832          call xmpi_sum(rhosu1_arr,spaceComm,ier)
1833          call xmpi_sum(rhosu2_arr,spaceComm,ier)
1834        end if
1835 
1836 !      Now symmetrize the magnetization at equivalent points
1837 !      Use Transpose[symrel]
1838        numpt=0
1839        do izone=1,izone_max
1840          rep=irrzon(izone,2,imagn)
1841          nup=nsym_used/rep
1842          do iup=1,nup
1843            ind=irrzon(iup+numpt,1,imagn)
1844 !          Get coordinates in the range [0,n-1]
1845            j=ind-1;j1=modulo(j,n1);j2=modulo(j/n1,n2);j3=j/(n1*n2)
1846 !TO BE UNCOMMENTED
1847 !          Get location of G vector (grid point) centered at 0 0 0
1848            l3=j3-(j3/id3)*n3
1849            l2=j2-(j2/id2)*n2
1850            l1=j1-(j1/id1)*n1
1851            if(fftn2_distrib(j2+1)==me_fft)  then ! this ind is to be treated by me_fft
1852              r2=ffti2_local(j2+1) - 1
1853              ind=n1*(nd2*j3+r2)+j1+1  ! this is ind in the current proc
1854              jsym=isymg(iup+numpt)
1855              if (jsym==0) then
1856                ABI_ERROR("jsym=0 in symrhg !")
1857              end if
1858              magxsu1=rhosu1_arr(3*izone-2);magxsu2=rhosu2_arr(3*izone-2)
1859              magysu1=rhosu1_arr(3*izone-1);magysu2=rhosu2_arr(3*izone-1)
1860              magzsu1=rhosu1_arr(3*izone  );magzsu2=rhosu2_arr(3*izone  )
1861 !            Working on this: the present coding will be detrimental for speed ! cos and sin are recomputed many times !
1862              tau1=tnons_used(1,jsym)
1863              tau2=tnons_used(2,jsym)
1864              tau3=tnons_used(3,jsym)
1865              if (abs(tau1)>tol12.or.abs(tau2)>tol12.or.abs(tau3)>tol12) then
1866 !              Compute exp(-2*Pi*I*G dot tau) using original G   (equivalent of phnons in the collinear case)
1867                arg=two_pi*(dble(l1)*tau1+dble(l2)*tau2+dble(l3)*tau3)
1868                phr=cos(arg)
1869                phi=-sin(arg)
1870              else
1871                phr=one
1872                phi=zero
1873              end if
1874              phr=phr*symafm_used(jsym)
1875              phi=phi*symafm_used(jsym)
1876 !TO BE COMMENTED
1877 !            phr=phnons(1,iup,imagn);if (rep==1) phr=phr*symafm_used(jsym) !if rep==2, symafm is already included in phnons
1878 !            phi=phnons(2,iup,imagn);if (rep==1) phi=phi*symafm_used(jsym) !(see irrzg.F90)
1879 !            The magnetization should transform as a vector in real space
1880 !            => symrel_cart  ?!?
1881              mxr=symrec_cart(1,1,jsym)*magxsu1+symrec_cart(2,1,jsym)*magysu1+symrec_cart(3,1,jsym)*magzsu1
1882              mxi=symrec_cart(1,1,jsym)*magxsu2+symrec_cart(2,1,jsym)*magysu2+symrec_cart(3,1,jsym)*magzsu2
1883              myr=symrec_cart(1,2,jsym)*magxsu1+symrec_cart(2,2,jsym)*magysu1+symrec_cart(3,2,jsym)*magzsu1
1884              myi=symrec_cart(1,2,jsym)*magxsu2+symrec_cart(2,2,jsym)*magysu2+symrec_cart(3,2,jsym)*magzsu2
1885              mzr=symrec_cart(1,3,jsym)*magxsu1+symrec_cart(2,3,jsym)*magysu1+symrec_cart(3,3,jsym)*magzsu1
1886              mzi=symrec_cart(1,3,jsym)*magxsu2+symrec_cart(2,3,jsym)*magysu2+symrec_cart(3,3,jsym)*magzsu2
1887 !            mxr=symrel_cart(1,1,jsym)*magxsu1+symrel_cart(1,2,jsym)*magysu1+symrel_cart(1,3,jsym)*magzsu1
1888 !            mxi=symrel_cart(1,1,jsym)*magxsu2+symrel_cart(1,2,jsym)*magysu2+symrel_cart(1,3,jsym)*magzsu2
1889 !            myr=symrel_cart(2,1,jsym)*magxsu1+symrel_cart(2,2,jsym)*magysu1+symrel_cart(2,3,jsym)*magzsu1
1890 !            myi=symrel_cart(2,1,jsym)*magxsu2+symrel_cart(2,2,jsym)*magysu2+symrel_cart(2,3,jsym)*magzsu2
1891 !            mzr=symrel_cart(3,1,jsym)*magxsu1+symrel_cart(3,2,jsym)*magysu1+symrel_cart(3,3,jsym)*magzsu1
1892 !            mzi=symrel_cart(3,1,jsym)*magxsu2+symrel_cart(3,2,jsym)*magysu2+symrel_cart(3,3,jsym)*magzsu2
1893              magngx(1,ind)=mxr*phr-mxi*phi
1894              magngx(2,ind)=mxi*phr+mxr*phi
1895              magngy(1,ind)=myr*phr-myi*phi
1896              magngy(2,ind)=myi*phr+myr*phi
1897              magngz(1,ind)=mzr*phr-mzi*phi
1898              magngz(2,ind)=mzi*phr+mzr*phi
1899            end if
1900          end do
1901          numpt=numpt+nup
1902        end do
1903        ABI_FREE(isymg)
1904        ABI_FREE(rhosu1_arr)
1905        ABI_FREE(rhosu2_arr)
1906        ABI_FREE(symrec_cart)
1907        ABI_FREE(symrel_cart)
1908        ABI_FREE(symafm_used)
1909        ABI_FREE(tnons_used)
1910 
1911      end if ! nspden==4
1912 
1913      call timab(17,2,tsec)
1914 
1915 !    Pull out full or spin up density, now symmetrized
1916      call fourdp(cplex,rhog,work,1,mpi_enreg,nfft,1,ngfft,0)
1917      rhor(:,ispden)=work(:)
1918      if (nspden==4) then
1919        call fourdp(cplex,magngx,work,1,mpi_enreg,nfft,1,ngfft,0)
1920        rhor(:,2)=work(:)
1921        call fourdp(cplex,magngy,work,1,mpi_enreg,nfft,1,ngfft,0)
1922        rhor(:,3)=work(:)
1923        call fourdp(cplex,magngz,work,1,mpi_enreg,nfft,1,ngfft,0)
1924        rhor(:,4)=work(:)
1925        ABI_FREE(magngx)
1926        ABI_FREE(magngy)
1927        ABI_FREE(magngz)
1928      end if
1929 
1930    end do ! ispden
1931 
1932  end if !  End on the condition nsym==1
1933 
1934  ABI_FREE(work)
1935 
1936  contains
1937 
1938    function map_symrhg(j1,n1)
1939 
1940    integer :: map_symrhg
1941    integer,intent(in) :: j1,n1
1942 !  Map into [0,n-1]
1943    map_symrhg=mod(n1+mod(j1,n1),n1)
1944  end function map_symrhg
1945 
1946 end subroutine symrhg