TABLE OF CONTENTS
- ABINIT/m_spacepar
- m_spacepar/hartre
- m_spacepar/hartrestr
- m_spacepar/irrzg
- m_spacepar/laplacian
- m_spacepar/meanvalue_g
- m_spacepar/redgr
- m_spacepar/rotate_rho
- m_spacepar/setsym
- m_spacepar/symrhg
ABINIT/m_spacepar [ Modules ]
NAME
m_spacepar
FUNCTION
Relatively Low-level procedures operating on arrays defined on the FFT box (G- or R- space) Unlike the procedures in m_cgtools, the routines declared in this module can use mpi_type.
COPYRIGHT
Copyright (C) 2008-2018 ABINIT group (XG, BA, MT, DRH, DCA, GMR, MJV) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt .
PARENTS
CHILDREN
SOURCE
22 #if defined HAVE_CONFIG_H 23 #include "config.h" 24 #endif 25 26 #include "abi_common.h" 27 28 module m_spacepar 29 30 use defs_basis 31 use m_abicore 32 use m_errors 33 use m_xmpi 34 use m_sort 35 36 use m_time, only : timab 37 use defs_abitypes, only : MPI_type 38 use m_symtk, only : mati3inv, chkgrp, symdet, symatm, matr3inv 39 use m_geometry, only : metric, symredcart 40 use m_mpinfo, only : ptabs_fourdp 41 use m_fft, only : zerosym, fourdp 42 43 implicit none 44 45 private
m_spacepar/hartre [ Functions ]
[ Top ] [ m_spacepar ] [ Functions ]
NAME
hartre
FUNCTION
Given rho(G), compute Hartree potential (=FFT of rho(G)/pi/(G+q)**2) When cplex=1, assume q=(0 0 0), and vhartr will be REAL When cplex=2, q must be taken into account, and vhartr will be COMPLEX
NOTES
*Modified code to avoid if statements inside loops to skip G=0. Replaced if statement on G^2>gsqcut to skip G s outside where rho(G) should be 0. Effect is negligible but gsqcut should be used to be strictly consistent with usage elsewhere in code. *The speed-up is provided by doing a few precomputations outside the inner loop. One variable size array is needed for this (gq).
INPUTS
cplex= if 1, vhartr is REAL, if 2, vhartr is COMPLEX gsqcut=cutoff value on G**2 for sphere inside fft box. (gsqcut=(boxcut**2)*ecut/(2.d0*(Pi**2)) izero=if 1, unbalanced components of Vhartree(g) are set to zero mpi_enreg=information about MPI parallelization nfft=(effective) number of FFT grid points (for this processor) ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft [qpt(3)=reduced coordinates for a wavevector to be combined with the G vectors (needed if cplex==2).] rhog(2,nfft)=electron density in G space rprimd(3,3)=dimensional primitive translations in real space (bohr) divgq0= [optional argument] value of the integration of the Coulomb singularity 4pi\int_BZ 1/q^2 dq
OUTPUT
vhartr(cplex*nfft)=Hartree potential in real space, either REAL or COMPLEX
PARENTS
dfpt_rhotov,dfptnl_loop,energy,fock_getghc,m_kxc,nonlinear,nres2vres odamix,prcref,prcref_PMA,respfn,rhotov,setup_positron,setvtr,tddft
CHILDREN
fourdp,metric,ptabs_fourdp,timab,zerosym
SOURCE
106 subroutine hartre(cplex,gsqcut,izero,mpi_enreg,nfft,ngfft,paral_kgb,rhog,rprimd,vhartr,& 107 & divgq0,qpt) ! Optional argument 108 109 110 !This section has been created automatically by the script Abilint (TD). 111 !Do not modify the following lines by hand. 112 #undef ABI_FUNC 113 #define ABI_FUNC 'hartre' 114 !End of the abilint section 115 116 implicit none 117 118 !Arguments ------------------------------------ 119 !scalars 120 integer,intent(in) :: cplex,izero,nfft,paral_kgb 121 real(dp),intent(in) :: gsqcut 122 type(MPI_type),intent(in) :: mpi_enreg 123 !arrays 124 integer,intent(in) :: ngfft(18) 125 real(dp),intent(in) :: rprimd(3,3),rhog(2,nfft) 126 real(dp),intent(in),optional :: divgq0 127 real(dp),intent(in),optional :: qpt(3) 128 real(dp),intent(out) :: vhartr(cplex*nfft) 129 130 !Local variables------------------------------- 131 !scalars 132 integer,parameter :: im=2,re=1 133 integer :: i1,i2,i23,i2_local,i3,id1,id2,id3 134 integer :: ig,ig1min,ig1,ig1max,ig2,ig2min,ig2max,ig3,ig3min,ig3max 135 integer :: ii,ii1,ing,n1,n2,n3,qeq0,qeq05,me_fft,nproc_fft 136 real(dp),parameter :: tolfix=1.000000001e0_dp 137 real(dp) :: cutoff,den,gqg2p3,gqgm12,gqgm13,gqgm23,gs,gs2,gs3,ucvol 138 character(len=500) :: message 139 !arrays 140 integer :: id(3) 141 integer, ABI_CONTIGUOUS pointer :: fftn2_distrib(:),ffti2_local(:) 142 integer, ABI_CONTIGUOUS pointer :: fftn3_distrib(:),ffti3_local(:) 143 real(dp) :: gmet(3,3),gprimd(3,3),qpt_(3),rmet(3,3),tsec(2) 144 real(dp),allocatable :: gq(:,:),work1(:,:) 145 146 ! ************************************************************************* 147 148 ! Keep track of total time spent in hartre 149 call timab(10,1,tsec) 150 151 ! Check that cplex has an allowed value 152 if(cplex/=1 .and. cplex/=2)then 153 write(message, '(a,i0,a,a)' )& 154 & 'From the calling routine, cplex=',cplex,ch10,& 155 & 'but the only value allowed are 1 and 2.' 156 MSG_BUG(message) 157 end if 158 159 call metric(gmet,gprimd,-1,rmet,rprimd,ucvol) 160 161 n1=ngfft(1); n2=ngfft(2); n3=ngfft(3) 162 nproc_fft = mpi_enreg%nproc_fft; me_fft = mpi_enreg%me_fft 163 164 ! Get the distrib associated with this fft_grid 165 call ptabs_fourdp(mpi_enreg,n2,n3,fftn2_distrib,ffti2_local,fftn3_distrib,ffti3_local) 166 167 ! Initialize a few quantities 168 cutoff=gsqcut*tolfix 169 if(present(qpt))then 170 qpt_=qpt 171 else 172 qpt_=zero 173 end if 174 qeq0=0 175 if(qpt_(1)**2+qpt_(2)**2+qpt_(3)**2<1.d-15) qeq0=1 176 qeq05=0 177 if (qeq0==0) then 178 if (abs(abs(qpt_(1))-half)<tol12.or.abs(abs(qpt_(2))-half)<tol12.or. & 179 & abs(abs(qpt_(3))-half)<tol12) qeq05=1 180 end if 181 182 ! If cplex=1 then qpt_ should be 0 0 0 183 if (cplex==1.and. qeq0/=1) then 184 write(message,'(a,3e12.4,a,a)')& 185 & 'cplex=1 but qpt=',qpt_,ch10,& 186 & 'qpt should be 0 0 0.' 187 MSG_BUG(message) 188 end if 189 190 ! If FFT parallelism then qpt should not be 1/2 191 if (nproc_fft>1.and.qeq05==1) then 192 write(message, '(a,3e12.4,a,a)' )& 193 & 'FFT parallelism selected but qpt',qpt_,ch10,& 194 & 'qpt(i) should not be 1/2...' 195 MSG_ERROR(message) 196 end if 197 198 ! In order to speed the routine, precompute the components of g+q 199 ! Also check if the booked space was large enough... 200 ABI_ALLOCATE(gq,(3,max(n1,n2,n3))) 201 do ii=1,3 202 id(ii)=ngfft(ii)/2+2 203 do ing=1,ngfft(ii) 204 ig=ing-(ing/id(ii))*ngfft(ii)-1 205 gq(ii,ing)=ig+qpt_(ii) 206 end do 207 end do 208 ig1max=-1;ig2max=-1;ig3max=-1 209 ig1min=n1;ig2min=n2;ig3min=n3 210 211 ABI_ALLOCATE(work1,(2,nfft)) 212 id1=n1/2+2;id2=n2/2+2;id3=n3/2+2 213 214 ! Triple loop on each dimension 215 do i3=1,n3 216 ig3=i3-(i3/id3)*n3-1 217 ! Precompute some products that do not depend on i2 and i1 218 gs3=gq(3,i3)*gq(3,i3)*gmet(3,3) 219 gqgm23=gq(3,i3)*gmet(2,3)*2 220 gqgm13=gq(3,i3)*gmet(1,3)*2 221 222 do i2=1,n2 223 ig2=i2-(i2/id2)*n2-1 224 if (fftn2_distrib(i2) == me_fft) then 225 gs2=gs3+ gq(2,i2)*(gq(2,i2)*gmet(2,2)+gqgm23) 226 gqgm12=gq(2,i2)*gmet(1,2)*2 227 gqg2p3=gqgm13+gqgm12 228 229 i2_local = ffti2_local(i2) 230 i23=n1*(i2_local-1 +(n2/nproc_fft)*(i3-1)) 231 ! Do the test that eliminates the Gamma point outside of the inner loop 232 ii1=1 233 if(i23==0 .and. qeq0==1 .and. ig2==0 .and. ig3==0)then 234 ii1=2 235 work1(re,1+i23)=zero 236 work1(im,1+i23)=zero 237 ! If the value of the integration of the Coulomb singularity 4pi\int_BZ 1/q^2 dq is given, use it 238 if (PRESENT(divgq0)) then 239 work1(re,1+i23)=rhog(re,1+i23)*divgq0*piinv 240 work1(im,1+i23)=rhog(im,1+i23)*divgq0*piinv 241 end if 242 end if 243 244 ! Final inner loop on the first dimension (note the lower limit) 245 do i1=ii1,n1 246 gs=gs2+ gq(1,i1)*(gq(1,i1)*gmet(1,1)+gqg2p3) 247 ii=i1+i23 248 249 if(gs<=cutoff)then 250 ! Identify min/max indexes (to cancel unbalanced contributions later) 251 ! Count (q+g)-vectors with similar norm 252 if ((qeq05==1).and.(izero==1)) then 253 ig1=i1-(i1/id1)*n1-1 254 ig1max=max(ig1max,ig1); ig1min=min(ig1min,ig1) 255 ig2max=max(ig2max,ig2); ig2min=min(ig2min,ig2) 256 ig3max=max(ig3max,ig3); ig3min=min(ig3min,ig3) 257 end if 258 259 den=piinv/gs 260 work1(re,ii)=rhog(re,ii)*den 261 work1(im,ii)=rhog(im,ii)*den 262 else 263 ! gs>cutoff 264 work1(re,ii)=zero 265 work1(im,ii)=zero 266 end if 267 268 end do ! End loop on i1 269 end if 270 end do ! End loop on i2 271 end do ! End loop on i3 272 273 ABI_DEALLOCATE(gq) 274 275 if (izero==1) then 276 ! Set contribution of unbalanced components to zero 277 278 if (qeq0==1) then !q=0 279 call zerosym(work1,2,n1,n2,n3,comm_fft=mpi_enreg%comm_fft,distribfft=mpi_enreg%distribfft) 280 281 else if (qeq05==1) then 282 !q=1/2; this doesn't work in parallel 283 ig1=-1;if (mod(n1,2)==0) ig1=1+n1/2 284 ig2=-1;if (mod(n2,2)==0) ig2=1+n2/2 285 ig3=-1;if (mod(n3,2)==0) ig3=1+n3/2 286 if (abs(abs(qpt_(1))-half)<tol12) then 287 if (abs(ig1min)<abs(ig1max)) ig1=abs(ig1max) 288 if (abs(ig1min)>abs(ig1max)) ig1=n1-abs(ig1min) 289 end if 290 if (abs(abs(qpt_(2))-half)<tol12) then 291 if (abs(ig2min)<abs(ig2max)) ig2=abs(ig2max) 292 if (abs(ig2min)>abs(ig2max)) ig2=n2-abs(ig2min) 293 end if 294 if (abs(abs(qpt_(3))-half)<tol12) then 295 if (abs(ig3min)<abs(ig3max)) ig3=abs(ig3max) 296 if (abs(ig3min)>abs(ig3max)) ig3=n3-abs(ig3min) 297 end if 298 call zerosym(work1,2,n1,n2,n3,ig1=ig1,ig2=ig2,ig3=ig3,& 299 & comm_fft=mpi_enreg%comm_fft,distribfft=mpi_enreg%distribfft) 300 end if 301 end if 302 303 ! Fourier Transform Vhartree. Vh in reciprocal space was stored in work1 304 call fourdp(cplex,work1,vhartr,1,mpi_enreg,nfft,ngfft,paral_kgb,0) 305 306 ABI_DEALLOCATE(work1) 307 308 call timab(10,2,tsec) 309 310 end subroutine hartre
m_spacepar/hartrestr [ Functions ]
[ Top ] [ m_spacepar ] [ Functions ]
NAME
hartrestr
FUNCTION
To be called for strain perturbation only Compute the inhomogenous terms generated by the strain derivative of Hartree potential due to the ground state charge rho(G) FFT of (rho(G)/pi)*[d(1/G**2)/d(strain) - delta(diagonal strain)*(1/G**2)]
NOTES
*based largely on hartre.f *Modified code to avoid if statements inside loops to skip G=0. Replaced if statement on G^2>gsqcut to skip G s outside where rho(G) should be 0. Effect is negligible but gsqcut should be used to be strictly consistent with usage elsewhere in code. *The speed-up is provided by doing a few precomputations outside the inner loop. One variable size array is needed for this (gq).
INPUTS
gsqcut=cutoff value on G**2 for sphere inside fft box. idir=direction of the current perturbation ipert=type of the perturbation mpi_enreg=information about MPI parallelization natom=number of atoms in cell. nfft=number of fft grid points (gsqcut=(boxcut**2)*ecut/(2._dp*(Pi**2)) ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft rhog(2,nfft)=array for Fourier transform of GS electron density rprimd(3,3)=dimensional primitive translations in real space (bohr)
OUTPUT
vhartr1(nfft)=Inhomogeneous term in strain-perturbation-induced Hartree potential in real space,
PARENTS
dfpt_nselt,dfpt_nstpaw,dfpt_rhotov
CHILDREN
fourdp,metric,ptabs_fourdp
SOURCE
911 subroutine hartrestr(gsqcut,idir,ipert,mpi_enreg,natom,nfft,ngfft,& 912 & paral_kgb,rhog,rprimd,vhartr1) 913 914 915 !This section has been created automatically by the script Abilint (TD). 916 !Do not modify the following lines by hand. 917 #undef ABI_FUNC 918 #define ABI_FUNC 'hartrestr' 919 !End of the abilint section 920 921 implicit none 922 923 !Arguments ------------------------------------ 924 !scalars 925 integer,intent(in) :: idir,ipert,natom,nfft,paral_kgb 926 real(dp),intent(in) :: gsqcut 927 type(MPI_type),intent(in) :: mpi_enreg 928 !arrays 929 integer,intent(in) :: ngfft(18) 930 real(dp),intent(in) :: rhog(2,nfft),rprimd(3,3) 931 real(dp),intent(out) :: vhartr1(nfft) 932 933 !Local variables------------------------------- 934 !scalars 935 integer,parameter :: im=2,re=1 936 integer :: i1,i2,i23,i3,id2,id3,ig,ig2,ig3,ii,ii1,ing,istr,ka,kb,n1,n2,n3 937 real(dp),parameter :: tolfix=1.000000001_dp 938 real(dp) :: cutoff,ddends,den,dgsds,gqg2p3,gqgm12,gqgm13,gqgm23,gs,gs2,gs3 939 real(dp) :: term,ucvol 940 character(len=500) :: message 941 !arrays 942 integer,save :: idx(12)=(/1,1,2,2,3,3,3,2,3,1,2,1/) 943 integer :: id(3) 944 integer, ABI_CONTIGUOUS pointer :: fftn2_distrib(:),ffti2_local(:) 945 integer, ABI_CONTIGUOUS pointer :: fftn3_distrib(:),ffti3_local(:) 946 real(dp) :: dgmetds(3,3),gmet(3,3),gprimd(3,3),gqr(3),rmet(3,3) 947 real(dp),allocatable :: gq(:,:),work1(:,:) 948 949 ! ************************************************************************* 950 951 if( .not. (ipert==natom+3 .or. ipert==natom+4))then 952 write(message, '(a,i0,a,a)' )& 953 & 'From the calling routine, ipert=',ipert,ch10,& 954 & 'so this routine for the strain perturbation should not be called.' 955 MSG_BUG(message) 956 end if 957 958 call metric(gmet,gprimd,-1,rmet,rprimd,ucvol) 959 960 n1=ngfft(1) ; n2=ngfft(2) ; n3=ngfft(3) 961 962 !Get the distrib associated with this fft_grid 963 call ptabs_fourdp(mpi_enreg,n2,n3,fftn2_distrib,ffti2_local,fftn3_distrib,ffti3_local) 964 965 !Initialize a few quantities 966 cutoff=gsqcut*tolfix 967 968 istr=idir + 3*(ipert-natom-3) 969 970 if(istr<1 .or. istr>6)then 971 write(message, '(a,i10,a,a,a)' )& 972 & 'Input dir gives istr=',istr,' not allowed.',ch10,& 973 & 'Possible values are 1,2,3,4,5,6 only.' 974 MSG_BUG(message) 975 end if 976 977 ka=idx(2*istr-1);kb=idx(2*istr) 978 do ii = 1,3 979 dgmetds(:,ii)=-(gprimd(ka,:)*gprimd(kb,ii)+gprimd(kb,:)*gprimd(ka,ii)) 980 end do 981 !For historical reasons: 982 dgmetds(:,:)=0.5_dp*dgmetds(:,:) 983 984 !In order to speed the routine, precompute the components of g+q 985 !Also check if the booked space was large enough... 986 ABI_ALLOCATE(gq,(3,max(n1,n2,n3))) 987 do ii=1,3 988 id(ii)=ngfft(ii)/2+2 989 do ing=1,ngfft(ii) 990 ig=ing-(ing/id(ii))*ngfft(ii)-1 991 gq(ii,ing)=ig 992 end do 993 end do 994 995 ABI_ALLOCATE(work1,(2,nfft)) 996 id2=n2/2+2 997 id3=n3/2+2 998 !Triple loop on each dimension 999 do i3=1,n3 1000 ig3=i3-(i3/id3)*n3-1 1001 ! Precompute some products that do not depend on i2 and i1 1002 gqr(3)=gq(3,i3) 1003 gs3=gq(3,i3)*gq(3,i3)*gmet(3,3) 1004 gqgm23=gq(3,i3)*gmet(2,3)*2 1005 gqgm13=gq(3,i3)*gmet(1,3)*2 1006 1007 do i2=1,n2 1008 if (fftn2_distrib(i2)==mpi_enreg%me_fft) then 1009 gqr(2)=gq(2,i2) 1010 gs2=gs3+ gq(2,i2)*(gq(2,i2)*gmet(2,2)+gqgm23) 1011 gqgm12=gq(2,i2)*gmet(1,2)*2 1012 gqg2p3=gqgm13+gqgm12 1013 ig2=i2-(i2/id2)*n2-1 1014 ! i23=n1*((i2-1)+n2*(i3-1)) 1015 i23=n1*((ffti2_local(i2)-1)+(n2/mpi_enreg%nproc_fft)*(i3-1)) 1016 ! Do the test that eliminates the Gamma point outside 1017 ! of the inner loop 1018 ii1=1 1019 if(i23==0 .and. ig2==0 .and. ig3==0)then 1020 ii1=2 1021 work1(re,1+i23)=0.0_dp 1022 work1(im,1+i23)=0.0_dp 1023 end if 1024 1025 ! Final inner loop on the first dimension 1026 ! (note the lower limit) 1027 do i1=ii1,n1 1028 gs=gs2+ gq(1,i1)*(gq(1,i1)*gmet(1,1)+gqg2p3) 1029 ii=i1+i23 1030 if(gs<=cutoff)then 1031 den=piinv/gs 1032 gqr(1)=gq(1,i1) 1033 dgsds=& 1034 & (gqr(1)*(dgmetds(1,1)*gqr(1)+dgmetds(1,2)*gqr(2)+dgmetds(1,3)*gqr(3))+ & 1035 & gqr(2)*(dgmetds(2,1)*gqr(1)+dgmetds(2,2)*gqr(2)+dgmetds(2,3)*gqr(3))+ & 1036 & gqr(3)*(dgmetds(3,1)*gqr(1)+dgmetds(3,2)*gqr(2)+dgmetds(3,3)*gqr(3)) ) 1037 ddends=-piinv*dgsds/gs**2 1038 if(istr<=3)then 1039 term=2.0_dp*ddends-den 1040 else 1041 term=2.0_dp*ddends 1042 end if 1043 work1(re,ii)=rhog(re,ii)*term 1044 work1(im,ii)=rhog(im,ii)*term 1045 else 1046 work1(re,ii)=0.0_dp 1047 work1(im,ii)=0.0_dp 1048 end if 1049 1050 end do ! End loop on i1 1051 end if 1052 end do ! End loop on i2 1053 end do ! End loop on i3 1054 1055 ABI_DEALLOCATE(gq) 1056 1057 !Fourier Transform Vhartree. 1058 !Vh in reciprocal space was stored in work1 1059 call fourdp(1,work1,vhartr1,1,mpi_enreg,nfft,ngfft,paral_kgb,0) 1060 1061 ABI_DEALLOCATE(work1) 1062 1063 end subroutine hartrestr
m_spacepar/irrzg [ Functions ]
[ Top ] [ m_spacepar ] [ Functions ]
NAME
irrzg
FUNCTION
Find the irreducible zone in reciprocal space under the symmetry group with real space rotations in symrel(3,3,nsym). The (integer) rotation matrices symrel(3,3,nsym) express the new real space positions (e.g. rotated atom positions) in REDUCED coordinates, i.e. in coordinates expressed as fractions of real space primitive translations (atomic coordinates xred). tnons(3,nsym) express the associated nonsymmorphic translations, again in reduced coordinates. Special data structure created in irrzon. First half holds mapping from irr zone to full zone; part of second half holds repetition number info. work1 is a work array to keep track of grid points found so far. In case nspden=2 and nsppol=1, one has to take care of antiferromagnetic operations. The subgroup of non-magnetic operations is used to generate irrzon(:,:,2) and phnons(:,:,2), while the full group is used to generate irrzon(:,:,1) and phnons(:,:,1)
NOTES
for reference in the near future (2018), some notes: this routine should be duplicated for
magnetizations in spinorial formalism. The only main difference will
be that the density is not simply transported to the image point under symrel
but is a vector which has to be transformed as well
$S \vec{m} (\vec{x}) = \sigma \vec{m} (S \vec{x} + \tau)$
$S \vec{m} (\vec{G}) = \sigma exp(+ 2 \pi i \vec{G} \vec{tau}) \vec{m}(S^{-1 T} \vec{G})$
S is a symop, sigma the AFM sign flip if any, tau the partial non symmorphic translation
x a position, m the magnetization 3 vector
The phase factor is the same as below for the density, but the collection of elements which
are equal is more complex: for a 3 or 6 axis the m vector could transform one component
to a linear combination of several others (I think). Things are not necessarily aligned
with the z axis.
INPUTS
nspden=number of spin-density components nsppol=1 for unpolarized, 2 for spin-polarized nsym=number of symmetry elements in group n1,n2,n3=box dimensions of real space grid (or fft box) symafm(nsym)=(anti)ferromagnetic part of symmetry operations symrel(3,3,nsym)=symmetry matrices in real space (integers) tnons(3,nsym)=reduced nonsymmorphic translations (symrel and tnons are in terms of real space primitive translations)
OUTPUT
irrzon(n1*n2*n3,2+(nspden/4),(nspden/nsppol)-3*(nspden/4))=integer array which contains the locations of related grid points and the repetition number for each symmetry class. phnons(2,n1*n2*n3,(nspden/nsppol)-3*(nspden/4))=phases associated with nonsymmorphic translations
PARENTS
m_ab7_kpoints,setsym,wfd_mkrho
CHILDREN
sort_int,wrtout
SOURCE
1575 subroutine irrzg(irrzon,nspden,nsppol,nsym,n1,n2,n3,phnons,symafm,symrel,tnons) 1576 1577 1578 !This section has been created automatically by the script Abilint (TD). 1579 !Do not modify the following lines by hand. 1580 #undef ABI_FUNC 1581 #define ABI_FUNC 'irrzg' 1582 !End of the abilint section 1583 1584 implicit none 1585 1586 !Arguments ------------------------------------ 1587 !scalars 1588 integer,intent(in) :: n1,n2,n3,nspden,nsppol,nsym 1589 !arrays 1590 integer,intent(in) :: symafm(nsym),symrel(3,3,nsym) 1591 integer,intent(out) :: irrzon(n1*n2*n3,2,(nspden/nsppol)-3*(nspden/4)) 1592 real(dp),intent(in) :: tnons(3,nsym) 1593 real(dp),intent(out) :: phnons(2,n1*n2*n3,(nspden/nsppol)-3*(nspden/4)) 1594 1595 !Local variables------------------------------- 1596 !scalars 1597 integer :: i1,i2,i3,id1,id2,id3,ifft,imagn,ind1,ind2,ipt,irep,isym,izone 1598 integer :: izonemax,j1,j2,j3,jj,k1,k2,k3,l1,l2,l3,nfftot,npt,nsym_used 1599 integer :: nzone,setzer,sppoldbl 1600 real(dp) :: arg,ph1i,ph1r,ph2i,ph2r,tau1,tau2,tau3 1601 logical,parameter :: afm_noncoll=.true. ! TRUE if antiferro symmetries are used in non-collinear magnetism 1602 character(len=500) :: message 1603 !arrays 1604 integer,allocatable :: class(:),iperm(:),symafm_used(:),symrel_used(:,:,:) 1605 integer,allocatable :: work1(:) 1606 real(dp),allocatable :: tnons_used(:,:),work2(:,:) 1607 1608 ! ************************************************************************* 1609 1610 ABI_ALLOCATE(class,(nsym)) 1611 ABI_ALLOCATE(iperm,(nsym)) 1612 ABI_ALLOCATE(work1,(n1*n2*n3)) 1613 ABI_ALLOCATE(work2,(2,n1*n2*n3)) 1614 1615 nfftot=n1*n2*n3 1616 1617 id1=n1/2+2 1618 id2=n2/2+2 1619 id3=n3/2+2 1620 1621 sppoldbl=nspden/nsppol;if (nspden==4) sppoldbl=1 1622 1623 do imagn=1,sppoldbl 1624 1625 ! Treat in a similar way the case of the full group and the non-magnetic subgroup 1626 nsym_used=0 1627 do isym=1,nsym 1628 if( (imagn==1 .and. sppoldbl==2) .or. symafm(isym)==1 .or. & 1629 & ((nspden==4).and.afm_noncoll) )then 1630 nsym_used=nsym_used+1 1631 end if 1632 end do 1633 1634 if(imagn==2 .and. nsym_used/=nsym/2)then 1635 write(message, '(a,a,a,a,a,i4,a,i0)' )& 1636 & ' The number of ferromagnetic symmetry operations must be',ch10,& 1637 & ' half the total number of operations, while it is observed that',ch10,& 1638 & ' nsym=',nsym,' and nsym_magn=',nsym_used 1639 MSG_BUG(message) 1640 end if 1641 1642 ABI_ALLOCATE(symafm_used,(nsym_used)) 1643 ABI_ALLOCATE(symrel_used,(3,3,nsym_used)) 1644 ABI_ALLOCATE(tnons_used,(3,nsym_used)) 1645 1646 nsym_used=0 1647 do isym=1,nsym 1648 if( (imagn==1 .and. sppoldbl==2) .or. symafm(isym)==1 .or. & 1649 & ((nspden==4).and.afm_noncoll) ) then 1650 nsym_used=nsym_used+1 1651 symrel_used(:,:,nsym_used)=symrel(:,:,isym) 1652 tnons_used(:,nsym_used)=tnons(:,isym) 1653 symafm_used(nsym_used)=symafm(isym) 1654 end if 1655 end do 1656 if ((nspden/=4).or.(.not.afm_noncoll)) symafm_used=1 1657 1658 1659 ! Zero out work array--later on, a zero entry will mean that 1660 ! a given grid point has not yet been assigned to an ibz point 1661 work1(1:nfftot)=0 1662 irrzon(:,2,imagn)=0 1663 1664 ! Initialize at G=0 (always in irreducible zone) 1665 nzone=1 1666 irrzon(1,1,imagn)=1 1667 irrzon(1,2,imagn)=nsym_used 1668 ! Set phase exp(2*Pi*I*G dot tnons) for G=0 1669 phnons(1,1,imagn)=one 1670 phnons(2,1,imagn)=zero 1671 npt=1 1672 ! setting work1(1)=1 indicates that first grid point (G=0) is 1673 ! in the iz (irreducible zone) 1674 work1(1)=1 1675 1676 ind1=0 1677 1678 ! Loop over reciprocal space grid points: 1679 do i3=1,n3 1680 do i2=1,n2 1681 do i1=1,n1 1682 1683 ind1=ind1+1 1684 1685 ! Check to see whether present grid point is equivalent to 1686 ! any previously identified ibz point--if not, a new ibz point 1687 ! has been found 1688 1689 if (work1(ind1)==0) then 1690 1691 ! A new class has been found. 1692 1693 ! Get location of G vector (grid point) centered at 0 0 0 1694 l3=i3-(i3/id3)*n3-1 1695 l2=i2-(i2/id2)*n2-1 1696 l1=i1-(i1/id1)*n1-1 1697 1698 do isym=1,nsym_used 1699 1700 ! Get rotated G vector Gj for each symmetry element 1701 ! -- here we use the TRANSPOSE of symrel_used; assuming symrel_used expresses 1702 ! the rotation in real space, the transpose is then appropriate 1703 ! for G space symmetrization (p. 1172d,e of notes, 2 June 1995). 1704 j1=symrel_used(1,1,isym)*l1+& 1705 & symrel_used(2,1,isym)*l2+symrel_used(3,1,isym)*l3 1706 j2=symrel_used(1,2,isym)*l1+& 1707 & symrel_used(2,2,isym)*l2+symrel_used(3,2,isym)*l3 1708 j3=symrel_used(1,3,isym)*l1+& 1709 & symrel_used(2,3,isym)*l2+symrel_used(3,3,isym)*l3 1710 1711 ! Map into [0,n-1] and then add 1 for array index in [1,n] 1712 k1=1+mod(n1+mod(j1,n1),n1) 1713 k2=1+mod(n2+mod(j2,n2),n2) 1714 k3=1+mod(n3+mod(j3,n3),n3) 1715 ! k1=1+map(j1,n1) 1716 ! k2=1+map(j2,n2) 1717 ! k3=1+map(j3,n3) 1718 1719 ! Get linear index of rotated point Gj 1720 ind2=k1+n1*((k2-1)+n2*(k3-1)) 1721 1722 ! Store info for new class: 1723 class(isym)=ind2 1724 iperm(isym)=isym 1725 1726 ! Setting work array element to 1 indicates grid point has been 1727 ! identified with iz point 1728 work1(ind2)=1 1729 1730 ! End of loop on isym 1731 end do 1732 1733 ! Sort integers into ascending order in each class 1734 ! (this lumps together G vectors with the same linear index, i.e. 1735 ! groups together symmetries which land on the same Gj) 1736 call sort_int(nsym_used,class,iperm) 1737 1738 ! Check repetition factor (how many identical copies of Gj occur 1739 ! from all symmetries applied to G) 1740 irep=0 1741 do isym=1,nsym_used 1742 if (class(isym)==class(1)) then 1743 irep=irep+1 1744 end if 1745 end do 1746 ipt=nsym_used/irep 1747 1748 ! Repetition factor must be divisor of nsym_used: 1749 if (nsym_used/=(ipt*irep)) then 1750 write(message, '(a,i5,a,i6,a,a,a,a,a,a)' )& 1751 & ' irep=',irep,' not a divisor of nsym_used=',nsym_used,ch10,& 1752 & ' This usually indicates that',& 1753 & ' the input symmetries do not form a group.',ch10,& 1754 & ' Action : check the input symmetries carefully do they',& 1755 & ' form a group ? If they do, there is a code bug.' 1756 MSG_ERROR(message) 1757 end if 1758 1759 ! Compute phases for any nonsymmorphic symmetries 1760 ! exp(-2*Pi*I*G dot tau(j)) for each symmetry j with 1761 ! (possibly zero) nonsymmorphic translation tau(j) 1762 do jj=1,nsym_used 1763 ! First get nonsymmorphic translation and see if nonzero 1764 ! (iperm grabs the symmetries in the new order after sorting) 1765 isym=iperm(jj) 1766 tau1=tnons_used(1,isym) 1767 tau2=tnons_used(2,isym) 1768 tau3=tnons_used(3,isym) 1769 if (abs(tau1)>tol12.or.abs(tau2)>tol12& 1770 & .or.abs(tau3)>tol12) then 1771 ! compute exp(-2*Pi*I*G dot tau) using original G 1772 arg=two_pi*(dble(l1)*tau1+dble(l2)*tau2+dble(l3)*tau3) 1773 work2(1,jj)=cos(arg) 1774 work2(2,jj)=-sin(arg) 1775 else 1776 work2(1,jj)=one 1777 work2(2,jj)=zero 1778 end if 1779 end do 1780 1781 ! All phases arising from symmetries which map to the same 1782 ! G vector must actually be the same because 1783 ! rho(Strans*G)=exp(2*Pi*I*(G) dot tau_S) rho(G) 1784 ! must be satisfied; if exp(2*Pi*I*(G) dot tau_S) can be different 1785 ! for two different symmetries S which both take G to the same St*G, 1786 ! then the related Fourier components rho(St*G) must VANISH. 1787 ! Thus: set "phase" to ZERO here in that case. 1788 ! The G mappings occur in sets of irep members; if irep=1 then 1789 ! all the G are unique. 1790 ! MT 090212: 1791 ! In the case of antiferromagn. symetries, one can have 1792 ! rho(Strans*G)= -exp(2*Pi*I*(G) dot tau_S) rho(G) 1793 ! (look at the minus !) 1794 ! A special treatment is then operated on phons. 1795 ! The later must be consistent with the use of phnons array 1796 ! in symrhg.F90 routine. 1797 ! XG 001108 : 1798 ! Note that there is a tolerance on the 1799 ! accuracy of tnons, especially when they are found from 1800 ! the symmetry finder (with xred that might be a bit inaccurate) 1801 if (irep > 1) then 1802 do jj=1,nsym_used,irep 1803 setzer=0 1804 ph1r=work2(1,jj);ph1i=work2(2,jj) 1805 do j1=jj,jj+irep-1 1806 ph2r=work2(1,j1);ph2i=work2(2,j1) 1807 if (((ph2r+ph1r)**2+(ph2i+ph1i)**2) <= tol14) then 1808 if (setzer/=1) setzer=-1 1809 else if (((ph2r-ph1r)**2+(ph2i-ph1i)**2) > tol14) then 1810 setzer=1 1811 end if 1812 end do 1813 ! Setzer= 0: phnons are all equal 1814 ! Setzer=-1: phnons are equal in absolute value 1815 ! Setzer= 1: some phnons are different 1816 if (setzer/=0) then 1817 if (setzer==-1) then 1818 if (afm_noncoll.and.nspden==4) then 1819 arg=symafm_used(iperm(jj)) 1820 if (all(symafm_used(iperm(jj:jj+irep-1))==arg)) then 1821 setzer=1 1822 else 1823 do j1=jj,jj+irep-1 1824 work2(:,j1)=work2(:,j1)*dble(symafm_used(iperm(j1))) 1825 end do 1826 end if 1827 else 1828 setzer=1 1829 end if 1830 end if 1831 if (setzer==1) work2(:,jj:jj+irep-1)=zero 1832 end if 1833 end do 1834 ! Compress data if irep>1: 1835 jj=0 1836 do isym=1,nsym_used,irep 1837 jj=jj+1 1838 class(jj)=class(isym) 1839 work2(1,jj)=work2(1,isym) 1840 work2(2,jj)=work2(2,isym) 1841 end do 1842 end if 1843 1844 ! Put new unique points into irrzon array: 1845 irrzon(1+npt:ipt+npt,1,imagn)=class(1:ipt) 1846 1847 ! Put repetition number into irrzon array: 1848 irrzon(1+nzone,2,imagn)=irep 1849 1850 ! DEBUG 1851 ! write(std_out,'(a,6i7)' )' irrzg : izone,i1,i2,i3,imagn,irrzon(859,2,1)=',& 1852 ! & 1+nzone,i1,i2,i3,imagn,irrzon(859,2,1) 1853 ! ENDDEBUG 1854 1855 ! Put phases (or 0) in phnons array: 1856 phnons(:,1+npt:ipt+npt,imagn)=work2(:,1:ipt) 1857 1858 ! Update number of points in irrzon array: 1859 ! (irep must divide evenly into nsym_used !) 1860 npt=npt+ipt 1861 1862 ! Update number of classes: 1863 nzone=nzone+1 1864 1865 end if 1866 ! 1867 ! End of loop on reciprocal space points, with indices i1, i2, i3 1868 end do 1869 end do 1870 end do 1871 1872 if (allocated(symafm_used)) then 1873 ABI_DEALLOCATE(symafm_used) 1874 end if 1875 if (allocated(symrel_used)) then 1876 ABI_DEALLOCATE(symrel_used) 1877 end if 1878 if (allocated(tnons_used)) then 1879 ABI_DEALLOCATE(tnons_used) 1880 end if 1881 1882 end do ! imagn 1883 1884 !Make sure number of real space points accounted for equals actual number of grid points 1885 if (npt/=n1*n2*n3) then 1886 write(message, '(a,a,a,a,i10,a,i10,a,a,a,a,a,a,a,a,a)' ) ch10,& 1887 & ' irrzg : ERROR -',ch10,& 1888 & ' npt=',npt,' and n1*n2*n3=',n1*n2*n3,' are not equal',ch10,& 1889 & ' This says that the total of all points in the irreducible',& 1890 & ' sector in real space',ch10,& 1891 & ' and all symmetrically equivalent',& 1892 & ' points, npt, does not equal the actual number',ch10,& 1893 & ' of real space grid points.' 1894 call wrtout(std_out,message,'COLL') 1895 write(message,'(3a)') & 1896 & ' This may mean that the input symmetries do not form a group',ch10,& 1897 & ' Action : check input symmetries carefully for errors.' 1898 MSG_ERROR(message) 1899 end if 1900 1901 !Perform some checks 1902 do imagn=1,sppoldbl 1903 1904 do ifft=1,nfftot 1905 if (irrzon(ifft,1,imagn)<1.or.irrzon(ifft,1,imagn)>nfftot) then 1906 write(message,'(a,4i0,a,a)')& 1907 & ' ifft,irrzon(ifft,1,imagn),nfftot,imagn=',ifft,irrzon(ifft,1,imagn),nfftot,imagn,ch10,& 1908 & ' =>irrzon goes outside acceptable bounds.' 1909 MSG_BUG(message) 1910 end if 1911 end do 1912 1913 izonemax=0 1914 do izone=1,nfftot 1915 ! Global bounds 1916 if (irrzon(izone,2,imagn)<0.or.irrzon(izone,2,imagn)>(nsym/imagn)) then 1917 write(message, '(a,5i7,a,a)' )& 1918 & ' izone,nzone,irrzon(izone,2,imagn),nsym,imagn =',izone,nzone,irrzon(izone,2,imagn),nsym,imagn,ch10,& 1919 & ' =>irrzon goes outside acceptable bounds.' 1920 MSG_BUG(message) 1921 end if 1922 ! Second index only goes up to nzone 1923 if(izonemax==0)then 1924 if (irrzon(izone,2,imagn)==0)izonemax=izone-1 1925 end if 1926 if(izonemax/=0)then 1927 if (irrzon(izone,2,imagn)/=0) then 1928 message = ' beyond izonemax, irrzon(izone,2,imagn) should be zero' 1929 MSG_BUG(message) 1930 end if 1931 end if 1932 end do 1933 1934 end do ! imagn 1935 1936 ABI_DEALLOCATE(class) 1937 ABI_DEALLOCATE(iperm) 1938 ABI_DEALLOCATE(work1) 1939 ABI_DEALLOCATE(work2) 1940 1941 end subroutine irrzg
m_spacepar/laplacian [ Functions ]
[ Top ] [ m_spacepar ] [ Functions ]
NAME
laplacian
FUNCTION
compute the laplacian of a function defined in real space the code is written in the way of /3xc/xcden.F90
INPUTS
gprimd(3,3)=dimensional reciprocal space primitive translations mpi_enreg=information about MPI parallelization nfft=number of points of the fft grid nfunc=number of functions on the grid for which the laplacian is to be calculated ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft paral_kgb=flag controlling (k,g,bands) parallelization (optional) rdfuncr(nfft,nfunc)=real(dp) discretized functions in real space rdfuncg_in TO BE DESCRIBED SB 090901 laplacerdfuncg_in TO BE DESCRIBED SB 090901 (optional) g2cart_in(nfft) = G**2 on the grid
OUTPUT
(optional) laplacerdfuncr = laplacian in real space of the functions in rdfuncr (optional) rdfuncg = real(dp) discretized functions in fourier space (optional) laplacerdfuncg = real(dp) discretized laplacian of the functions in fourier space (optional) g2cart_out(nfft) = G**2 on the grid rdfuncg_out TO BE DESCRIBED SB 090901 laplacerdfuncg_out TO BE DESCRIBED SB 090901
PARENTS
frskerker1,frskerker2,moddiel_csrb,prcrskerker1,prcrskerker2
CHILDREN
fourdp,ptabs_fourdp
SOURCE
565 subroutine laplacian(gprimd,mpi_enreg,nfft,nfunc,ngfft,paral_kgb,rdfuncr,& 566 & laplacerdfuncr,rdfuncg_out,laplacerdfuncg_out,g2cart_out,rdfuncg_in,g2cart_in) 567 568 569 !This section has been created automatically by the script Abilint (TD). 570 !Do not modify the following lines by hand. 571 #undef ABI_FUNC 572 #define ABI_FUNC 'laplacian' 573 !End of the abilint section 574 575 implicit none 576 577 !Arguments ------------------------------------ 578 !scalars 579 integer,intent(in) :: nfft,nfunc,paral_kgb 580 type(MPI_type),intent(in) :: mpi_enreg 581 !arrays 582 integer,intent(in) :: ngfft(18) 583 real(dp),intent(in) :: gprimd(3,3) 584 real(dp),intent(inout),optional :: laplacerdfuncr(nfft,nfunc) 585 real(dp),intent(inout),optional,target :: rdfuncr(nfft,nfunc) 586 real(dp),intent(in),optional,target :: g2cart_in(nfft) !vz_i 587 real(dp),intent(out),optional,target :: g2cart_out(nfft) !vz_i 588 real(dp),intent(out),optional,target :: laplacerdfuncg_out(2,nfft,nfunc) 589 real(dp),intent(in),optional,target :: rdfuncg_in(2,nfft,nfunc) !vz_i 590 real(dp),intent(out),optional,target :: rdfuncg_out(2,nfft,nfunc) 591 592 !Local variables------------------------------- 593 !scalars 594 integer :: count,i1,i2,i3,id1,id2,id3,ifft,ifunc,ig1,ig2,ig3,ii1,n1,n2 595 integer :: n3 596 real(dp) :: b11,b12,b13,b21,b22,b23,b31,b32,b33 597 !arrays 598 integer, ABI_CONTIGUOUS pointer :: fftn2_distrib(:),ffti2_local(:) 599 integer, ABI_CONTIGUOUS pointer :: fftn3_distrib(:),ffti3_local(:) 600 real(dp),ABI_CONTIGUOUS pointer :: g2cart(:),laplacerdfuncg(:,:,:),rdfuncg(:,:,:) 601 602 ! ************************************************************************* 603 604 !Keep local copy of fft dimensions 605 n1=ngfft(1) 606 n2=ngfft(2) 607 n3=ngfft(3) 608 609 if(present(laplacerdfuncg_out)) then 610 laplacerdfuncg => laplacerdfuncg_out 611 else 612 ABI_ALLOCATE(laplacerdfuncg,(2,nfft,nfunc)) 613 end if 614 615 ! Get the distrib associated with this fft_grid 616 call ptabs_fourdp(mpi_enreg,n2,n3,fftn2_distrib,ffti2_local,fftn3_distrib,ffti3_local) 617 618 !change the real density rdfuncr on real space on the real density 619 !rdfuncg in reciprocal space 620 if(.not.present(rdfuncg_in)) then 621 if(present(rdfuncg_out)) then 622 rdfuncg => rdfuncg_out 623 else 624 ABI_ALLOCATE(rdfuncg,(2,nfft,nfunc)) 625 end if 626 if(present(rdfuncr)) then 627 do ifunc=1,nfunc 628 call fourdp(1,rdfuncg(:,:,ifunc),rdfuncr(:,ifunc),-1,mpi_enreg,nfft,ngfft,paral_kgb,0) 629 end do 630 end if 631 else 632 rdfuncg => rdfuncg_in 633 end if 634 635 !apply the laplacian on laplacerdfuncr 636 !code from /3xc/xcden.F90 637 !see also src/5common/hatre.F90 and src/5common/moddiel.F90 638 !Keep local copy of fft dimensions 639 !Initialize computation of G^2 in cartesian coordinates 640 if(.not.present(g2cart_in)) then 641 if(present(g2cart_out)) then 642 g2cart => g2cart_out 643 else 644 ABI_ALLOCATE(g2cart,(nfft)) 645 end if 646 id1=int(n1/2)+2 647 id2=int(n2/2)+2 648 id3=int(n3/2)+2 649 count=0 650 do i3=1,n3 651 ifft=(i3-1)*n1*(n2/mpi_enreg%nproc_fft) 652 ig3=i3-int(i3/id3)*n3-1 653 do i2=1,n2 654 if (fftn2_distrib(i2)==mpi_enreg%me_fft) then 655 ig2=i2-int(i2/id2)*n2-1 656 657 ii1=1 658 do i1=ii1,n1 659 ig1=i1-int(i1/id1)*n1-1 660 ifft=ifft+1 661 662 b11=gprimd(1,1)*real(ig1,dp) 663 b21=gprimd(2,1)*real(ig1,dp) 664 b31=gprimd(3,1)*real(ig1,dp) 665 b12=gprimd(1,2)*real(ig2,dp) 666 b22=gprimd(2,2)*real(ig2,dp) 667 b32=gprimd(3,2)*real(ig2,dp) 668 b13=gprimd(1,3)*real(ig3,dp) 669 b23=gprimd(2,3)*real(ig3,dp) 670 b33=gprimd(3,3)*real(ig3,dp) 671 672 g2cart(ifft)=( & 673 & (b11+b12+b13)**2& 674 & +(b21+b22+b23)**2& 675 & +(b31+b32+b33)**2& 676 & ) 677 do ifunc=1,nfunc 678 ! compute the laplacian in Fourier space that is * (i x 2pi x G)**2 679 laplacerdfuncg(1,ifft,ifunc) = -rdfuncg(1,ifft,ifunc)*g2cart(ifft)*two_pi*two_pi 680 laplacerdfuncg(2,ifft,ifunc) = -rdfuncg(2,ifft,ifunc)*g2cart(ifft)*two_pi*two_pi 681 end do 682 end do 683 end if 684 end do 685 end do 686 if(.not.present(g2cart_out)) then 687 ABI_DEALLOCATE(g2cart) 688 end if 689 else 690 g2cart => g2cart_in 691 do ifunc=1,nfunc 692 do ifft=1,nfft 693 ! compute the laplacian in Fourier space that is * (i x 2pi x G)**2 694 laplacerdfuncg(1,ifft,ifunc) = -rdfuncg(1,ifft,ifunc)*g2cart(ifft)*two_pi*two_pi 695 laplacerdfuncg(2,ifft,ifunc) = -rdfuncg(2,ifft,ifunc)*g2cart(ifft)*two_pi*two_pi 696 end do 697 end do 698 end if 699 700 !get the result back into real space 701 if(present(laplacerdfuncr)) then 702 do ifunc=1,nfunc 703 call fourdp(1,laplacerdfuncg(:,:,ifunc),laplacerdfuncr(:,ifunc),1,mpi_enreg,nfft,ngfft,paral_kgb,0) 704 end do 705 end if 706 707 !deallocate pointers 708 if((.not.present(rdfuncg_in)).and.(.not.present(rdfuncg_in))) then 709 ABI_DEALLOCATE(rdfuncg) 710 end if 711 if(.not.present(laplacerdfuncg_out)) then 712 ABI_DEALLOCATE(laplacerdfuncg) 713 end if 714 715 end subroutine laplacian
m_spacepar/meanvalue_g [ Functions ]
[ Top ] [ m_spacepar ] [ Functions ]
NAME
meanvalue_g
FUNCTION
Compute the mean value of one wavefunction, in reciprocal space, for an operator that is real, diagonal in G-space: <wf|op|wf> For the time being, only spin-independent operators are treated.
INPUTS
diag(npw)=diagonal operator (real, spin-independent!)
filter= if 1, need to filter on the value of diag, that must be less than huge(0.0d0)*1.d-11
otherwise, should be 0
istwf_k=storage mode of the vectors
npw=number of planewaves of the vector
nspinor=number of spinor components
vect(2,npw*nspinor)=vector
vect1(2,npw*nspinor*use_ndo)=vector1 (=vector in most of the cases)
use_ndo = says if vect=/vect1
OUTPUT
ar=mean value
PARENTS
dfpt_vtowfk,energy,forstrnps,vtowfk
CHILDREN
xmpi_sum
SOURCE
344 subroutine meanvalue_g(ar,diag,filter,istwf_k,mpi_enreg,npw,nspinor,vect,vect1,use_ndo,ar_im) 345 346 347 !This section has been created automatically by the script Abilint (TD). 348 !Do not modify the following lines by hand. 349 #undef ABI_FUNC 350 #define ABI_FUNC 'meanvalue_g' 351 !End of the abilint section 352 353 implicit none 354 355 !Arguments ------------------------------------ 356 !scalars 357 integer,intent(in) :: filter,istwf_k,npw,nspinor,use_ndo 358 real(dp),intent(out) :: ar 359 real(dp),intent(out),optional :: ar_im 360 type(MPI_type),intent(in) :: mpi_enreg 361 !arrays 362 real(dp),intent(in) :: diag(npw),vect(2,npw*nspinor) 363 real(dp),intent(in) :: vect1(2,npw*nspinor) 364 365 !Local variables------------------------------- 366 !scalars 367 integer :: i1,ierr,ipw,jpw,me_g0 368 character(len=500) :: message 369 !arrays 370 371 ! ************************************************************************* 372 me_g0 = mpi_enreg%me_g0 373 374 DBG_CHECK(ANY(filter==(/0,1/)),"Wrong filter") 375 DBG_CHECK(ANY(nspinor==(/1,2/)),"Wrong nspinor") 376 DBG_CHECK(ANY(istwf_k==(/(ipw,ipw=1,9)/)),"Wrong istwf_k") 377 378 if(nspinor==2 .and. istwf_k/=1)then 379 write(message,'(a,a,a,i6,a,i6)')& 380 & 'When istwf_k/=1, nspinor must be 1,',ch10,& 381 & 'however, nspinor=',nspinor,', and istwf_k=',istwf_k 382 MSG_BUG(message) 383 end if 384 385 if(use_ndo==1 .and. (istwf_k==2 .and.me_g0==1)) then 386 MSG_BUG('use_ndo==1, not tested') 387 end if 388 389 ar=zero 390 if(present(ar_im)) ar_im=zero 391 392 !Normal storage mode 393 if(istwf_k==1)then 394 395 ! No filter 396 if(filter==0)then 397 !$OMP PARALLEL DO REDUCTION(+:ar) 398 do ipw=1,npw 399 ar=ar+diag(ipw)*(vect(1,ipw)*vect1(1,ipw)+vect(2,ipw)*vect1(2,ipw)) 400 end do 401 if(nspinor==2)then 402 !$OMP PARALLEL DO REDUCTION(+:ar) PRIVATE(jpw) 403 do ipw=1+npw,2*npw 404 jpw=ipw-npw 405 ar=ar+diag(jpw)*(vect(1,ipw)*vect1(1,ipw)+vect(2,ipw)*vect1(2,ipw)) 406 end do 407 end if 408 if(use_ndo==1 .and. nspinor==2)then 409 !$OMP PARALLEL DO REDUCTION(+:ar_im) 410 do ipw=1,npw 411 ar_im=ar_im+diag(ipw)*(vect1(1,ipw)*vect(2,ipw)-vect1(2,ipw)*vect(1,ipw)) 412 end do 413 !$OMP PARALLEL DO REDUCTION(+:ar_im) PRIVATE(jpw) 414 do ipw=1+npw,2*npw 415 jpw=ipw-npw 416 ar_im=ar_im+diag(jpw)*(vect1(1,ipw)*vect(2,ipw)-vect1(2,ipw)*vect(1,ipw)) 417 end do 418 end if 419 420 ! !$OMP PARALLEL DO REDUCTION(+:ar,ar_im) 421 ! do ipw=1,npw 422 ! ar=ar+diag(ipw)*(vect(1,ipw)*vect1(1,ipw)+vect(2,ipw)*vect1(2,ipw)) 423 ! if(use_ndo==1.and.nspinor==2) ar_im=ar_im+diag(ipw)*(vect1(1,ipw)*vect(2,ipw)-vect1(2,ipw)*vect(1,ipw)) 424 ! end do 425 ! if(nspinor==2)then 426 ! !$OMP PARALLEL DO PRIVATE(ipw) REDUCTION(+:ar,ar_im) 427 ! do ipw=1+npw,2*npw 428 ! ar=ar+diag(ipw-npw)*(vect(1,ipw)*vect1(1,ipw)+vect(2,ipw)*vect1(2,ipw)) 429 ! if(use_ndo==1.and.nspinor==2) ar_im=ar_im+diag(ipw-npw)*(vect1(1,ipw)*vect(2,ipw)-vect1(2,ipw)*vect(1,ipw)) 430 ! end do 431 ! end if 432 else ! will filter 433 434 !$OMP PARALLEL DO REDUCTION(+:ar) 435 do ipw=1,npw 436 if(diag(ipw)<huge(0.0d0)*1.d-11)then 437 ar=ar+diag(ipw)*(vect(1,ipw)*vect1(1,ipw)+vect(2,ipw)*vect1(2,ipw)) 438 end if 439 end do 440 if(nspinor==2)then 441 !$OMP PARALLEL DO REDUCTION(+:ar) PRIVATE(jpw) 442 do ipw=1+npw,2*npw 443 jpw=ipw-npw 444 if(diag(jpw)<huge(0.0d0)*1.d-11)then 445 ar=ar+diag(jpw)*(vect(1,ipw)*vect1(1,ipw)+vect(2,ipw)*vect1(2,ipw)) 446 end if 447 end do 448 end if 449 if(use_ndo==1 .and. nspinor==2)then 450 !$OMP PARALLEL DO REDUCTION(+:ar_im) 451 do ipw=1,npw 452 if(diag(ipw)<huge(0.0d0)*1.d-11)then 453 ar_im=ar_im+diag(ipw)*(vect1(1,ipw)*vect(2,ipw)-vect1(2,ipw)*vect(1,ipw)) 454 end if 455 end do 456 !$OMP PARALLEL DO REDUCTION(+:ar_im) PRIVATE(jpw) 457 do ipw=1+npw,2*npw 458 jpw=ipw-npw 459 if(diag(jpw)<huge(0.0d0)*1.d-11)then 460 ar_im=ar_im+diag(jpw)*(vect1(1,ipw)*vect(2,ipw)-vect1(2,ipw)*vect(1,ipw)) 461 end if 462 end do 463 end if 464 465 466 ! !$OMP PARALLEL DO PRIVATE(ipw) REDUCTION(+:ar,ar_im) 467 ! do ipw=1,npw 468 ! if(diag(ipw)<huge(0.0d0)*1.d-11)then 469 ! ar=ar+diag(ipw)*(vect(1,ipw)*vect1(1,ipw)+vect(2,ipw)*vect1(2,ipw)) 470 ! if(use_ndo==1.and.nspinor==2) ar_im=ar_im+diag(ipw)*(vect1(1,ipw)*vect(2,ipw)-vect1(2,ipw)*vect(1,ipw)) 471 ! end if 472 ! end do 473 ! if(nspinor==2)then 474 ! !$OMP PARALLEL DO PRIVATE(ipw) REDUCTION(+:ar,ar_im) 475 ! do ipw=1+npw,2*npw 476 ! if(diag(ipw-npw)<huge(0.0d0)*1.d-11)then 477 ! ar=ar+diag(ipw-npw)*(vect(1,ipw)*vect1(1,ipw)+vect(2,ipw)*vect1(2,ipw)) 478 ! if(use_ndo==1.and.nspinor==2) ar_im=ar_im+diag(ipw-npw)*(vect1(1,ipw)*vect(2,ipw)-vect1(2,ipw)*vect(1,ipw)) 479 ! end if 480 ! end do 481 ! end if ! nspinor==2 482 483 end if ! filter==0 484 485 else if(istwf_k>=2)then 486 487 if(filter==0)then 488 i1=1 489 if(istwf_k==2 .and. me_g0==1)then ! MPIWF need to know which proc has G=0 490 ar=half*diag(1)*vect(1,1)*vect1(1,1) ; i1=2 491 end if 492 493 !$OMP PARALLEL DO REDUCTION(+:ar) 494 do ipw=i1,npw 495 ar=ar+diag(ipw)*(vect(1,ipw)*vect1(1,ipw)+vect(2,ipw)*vect1(2,ipw)) 496 end do 497 498 else ! filter/=0 499 i1=1 500 if(istwf_k==2 .and. me_g0==1)then 501 if(diag(1)<huge(0.0d0)*1.d-11)then 502 ar=half*diag(1)*vect(1,1)*vect1(1,1) ; i1=2 503 end if 504 end if 505 506 !$OMP PARALLEL DO REDUCTION(+:ar) 507 do ipw=i1,npw 508 if(diag(ipw)<huge(0.0d0)*1.d-11)then 509 ar=ar+diag(ipw)*(vect(1,ipw)*vect1(1,ipw)+vect(2,ipw)*vect1(2,ipw)) 510 end if 511 end do 512 end if ! filter==0 513 514 ar=two*ar 515 516 end if ! istwf_k 517 518 !MPIWF need to make reduction on ar and ai . 519 if(mpi_enreg%paral_kgb==1)then 520 call xmpi_sum(ar,mpi_enreg%comm_bandspinorfft ,ierr) 521 if(present(ar_im))then 522 call xmpi_sum(ar_im,mpi_enreg%comm_bandspinorfft,ierr) 523 end if 524 end if 525 526 end subroutine meanvalue_g
m_spacepar/redgr [ Functions ]
[ Top ] [ m_spacepar ] [ Functions ]
NAME
redgr
FUNCTION
Compute reduced gradients of a real function on the usual unshifted fft grid. The gradient directions are the along the primitive reciprocal lattice vectors. The input function is intended to be a single spin component of the valence charge density, the valence + core charge densities or the first-order core charge density for use in frozen wf elastic tensor calculations within the GGA.
NOTES
Closely linked to xcden, but limited to Q=0, real charge densities, and unshifted grids.
INPUTS
mpi_enreg=information about MPI parallelization nfft=(effective) number of FFT grid points (for this processor) ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft frin(nfft)=real space input function
OUTPUT
frredgr(nfft,3)= reduced gradient of input function (same units as frin)
PARENTS
dfpt_eltfrxc
CHILDREN
fourdp
SOURCE
752 subroutine redgr (frin,frredgr,mpi_enreg,nfft,ngfft,paral_kgb) 753 754 755 !This section has been created automatically by the script Abilint (TD). 756 !Do not modify the following lines by hand. 757 #undef ABI_FUNC 758 #define ABI_FUNC 'redgr' 759 !End of the abilint section 760 761 implicit none 762 763 !Arguments ------------------------------------ 764 !scalars 765 integer,intent(in) :: nfft,paral_kgb 766 type(MPI_type),intent(in) :: mpi_enreg 767 !arrays 768 integer,intent(in) :: ngfft(18) 769 real(dp),intent(in) :: frin(nfft) 770 real(dp),intent(out) :: frredgr(nfft,3) 771 772 !Local variables------------------------------- 773 !scalars 774 integer :: cplex_tmp,i1,i2,i3,id,idir,ifft,ig,ii,ing,n1,n2,n3 775 !arrays 776 real(dp),allocatable :: gg(:,:),wkcmpx(:,:),work(:),workgr(:,:) 777 778 ! ************************************************************************* 779 780 !Only real arrays are treated 781 cplex_tmp=1 782 783 !Keep local copy of fft dimensions 784 n1=ngfft(1) ; n2=ngfft(2) ; n3=ngfft(3) 785 786 !In order to speed the routine, precompute the components of g, including 2pi factor 787 ABI_ALLOCATE(gg,(max(n1,n2,n3),3)) 788 do ii=1,3 789 id=ngfft(ii)/2+2 790 do ing=1,ngfft(ii) 791 ig=ing-(ing/id)*ngfft(ii)-1 792 gg(ing,ii)=two_pi*ig 793 end do 794 ! Note that the G <-> -G symmetry must be maintained 795 if(mod(ngfft(ii),2)==0)gg(ngfft(ii)/2+1,ii)=zero 796 end do 797 798 ABI_ALLOCATE(wkcmpx,(2,nfft)) 799 ABI_ALLOCATE(work,(nfft)) 800 ABI_ALLOCATE(workgr,(2,nfft)) 801 802 !Obtain rho(G) in wkcmpx from input rho(r) 803 work(:)=frin(:) 804 805 call fourdp(cplex_tmp,wkcmpx,work,-1,mpi_enreg,nfft,ngfft,paral_kgb,0) 806 807 !Gradient calculation for three reduced components in turn. 808 !Code duplicated to remove logic from loops. 809 do idir=1,3 810 if(idir==1) then 811 !$OMP PARALLEL DO PRIVATE(ifft) 812 do i3=1,n3 813 ifft=(i3-1)*n1*n2 814 do i2=1,n2 815 do i1=1,n1 816 ifft=ifft+1 817 ! Multiply by i 2pi G(idir) 818 workgr(2,ifft)= gg(i1,idir)*wkcmpx(1,ifft) 819 workgr(1,ifft)=-gg(i1,idir)*wkcmpx(2,ifft) 820 end do 821 end do 822 end do 823 else if(idir==2) then 824 !$OMP PARALLEL DO PRIVATE(ifft) 825 do i3=1,n3 826 ifft=(i3-1)*n1*n2 827 do i2=1,n2 828 do i1=1,n1 829 ifft=ifft+1 830 ! Multiply by i 2pi G(idir) 831 workgr(2,ifft)= gg(i2,idir)*wkcmpx(1,ifft) 832 workgr(1,ifft)=-gg(i2,idir)*wkcmpx(2,ifft) 833 end do 834 end do 835 end do 836 else 837 !$OMP PARALLEL DO PRIVATE(ifft) 838 do i3=1,n3 839 ifft=(i3-1)*n1*n2 840 do i2=1,n2 841 do i1=1,n1 842 ifft=ifft+1 843 ! Multiply by i 2pi G(idir) 844 workgr(2,ifft)= gg(i3,idir)*wkcmpx(1,ifft) 845 workgr(1,ifft)=-gg(i3,idir)*wkcmpx(2,ifft) 846 end do 847 end do 848 end do 849 end if !idir 850 851 call fourdp(cplex_tmp,workgr,work,1,mpi_enreg,nfft,ngfft,paral_kgb,0) 852 853 !$OMP PARALLEL DO 854 do ifft=1,nfft 855 frredgr(ifft,idir)=work(ifft) 856 end do 857 858 end do !idir 859 860 ABI_DEALLOCATE(gg) 861 ABI_DEALLOCATE(wkcmpx) 862 ABI_DEALLOCATE(work) 863 ABI_DEALLOCATE(workgr) 864 865 end subroutine redgr
m_spacepar/rotate_rho [ Functions ]
[ Top ] [ m_spacepar ] [ Functions ]
NAME
rotate_rho
FUNCTION
rotate density in real and reciprocal space
INPUTS
cplex: if 1, real space 1-order functions on FFT grid are REAL, if 2, COMPLEX mpi_enreg=informations about MPI parallelization nfft=(effective) number of FFT grid points (for this proc) for the "fine" grid (see NOTES in respfn.F90) ngfft=array of dimensions for different FFT grids nspden=number of spin-density components rhor1(cplex*nfft,nspden)=array for Fourier transform of RF electron density === if psps%usepaw==1 TODO: extend to PAW pawrhoij1(natom) <type(pawrhoij_type)>= 1st-order paw rhoij occupancies and related data symrel1=single symmetry operation in real space to apply to rho tnon = eventual translation associated to symrel1
OUTPUT
rhog1_eq(2,nfft)= symmetric density in reciprocal space for equivalent perturbation rhor1_eq(cplex*nfft,nspden) = symmetric density in real space for equivalent perturbation
PARENTS
dfpt_looppert
CHILDREN
fourdp
TODO
This routine is deprecated. Now the symmetrization of the **potentials** is done in the m_dvdb Moreover one should take linear combinations of symmetrical components
SOURCE
1979 subroutine rotate_rho(cplex, itirev, mpi_enreg, nfft, ngfft, nspden, & 1980 & rhor1, rhog1_eq, rhor1_eq, symrel1, tnon) 1981 1982 1983 !This section has been created automatically by the script Abilint (TD). 1984 !Do not modify the following lines by hand. 1985 #undef ABI_FUNC 1986 #define ABI_FUNC 'rotate_rho' 1987 !End of the abilint section 1988 1989 implicit none 1990 1991 !args 1992 integer,intent(in) :: cplex, nfft, nspden, itirev 1993 integer,intent(in) :: ngfft(18) 1994 1995 integer, intent(in) :: symrel1(3,3) 1996 real(dp),intent(in) :: tnon(3) 1997 real(dp),intent(inout) :: rhor1(cplex*nfft,nspden) 1998 1999 real(dp),intent(out) :: rhog1_eq(2,nfft) 2000 real(dp),intent(out) :: rhor1_eq(cplex*nfft,nspden) 2001 2002 type(MPI_type),intent(in) :: mpi_enreg 2003 2004 ! local vars 2005 integer :: id1, id2, id3 2006 integer :: n1, n2, n3, nd2 2007 integer :: l1, l2, l3 2008 integer :: i1, i2, i3, ind1, ind2 2009 integer :: j1, j2, j3 2010 integer :: k1, k2, k3 2011 integer :: nproc_fft, ispden, me_fft 2012 real(dp) :: arg 2013 logical :: t_tnon_nonzero 2014 2015 real(dp) :: phnon1(2) 2016 real(dp), allocatable :: workg(:,:), workg_eq(:,:) 2017 character(len=500) :: message 2018 2019 ! ************************************************************************* 2020 2021 n1=ngfft(1);n2=ngfft(2);n3=ngfft(3);nproc_fft=ngfft(10);me_fft=ngfft(11);nd2=n2/nproc_fft 2022 2023 id1=n1/2+2 2024 id2=n2/2+2 2025 id3=n3/2+2 2026 2027 rhog1_eq = zero 2028 rhor1_eq = zero 2029 2030 if (itirev == 2) then 2031 write (message,'(3a,9I4,1a)') 'using time reversal. ',ch10,'Symrel1 = ', symrel1, ch10 2032 else 2033 write (message,'(3a,9I4,1a)') 'no time reversal. ',ch10,'Symrel1 = ', symrel1, ch10 2034 end if 2035 !call wrtout(std_out, message, 'COLL') 2036 2037 t_tnon_nonzero = (any(abs(tnon) > tol12)) 2038 2039 ! eventually, for FFT parallelization 2040 ! call ptabs_fourdp(mpi_enreg,n2,n3,fftn2_distrib,ffti2_local,fftn3_distrib,ffti3_local) 2041 2042 ABI_ALLOCATE(workg,(2,nfft)) 2043 ABI_ALLOCATE(workg_eq,(2,nfft)) 2044 do ispden = 1, nspden 2045 2046 ! fft input rhor1 to reciprocal space: uses work* as a buffer 2047 call fourdp(cplex,workg,rhor1(:,ispden),-1,mpi_enreg,nfft,ngfft,mpi_enreg%paral_kgb,0) 2048 2049 ! below taken from irrzg and setsym 2050 ! Loop over reciprocal space grid points: 2051 ! loop over local points in workg, and get back transform from rhog1, 2052 ! which is presumed complete on each proc 2053 ind1=0 2054 do i3=1,n3 2055 do i2=1,n2 2056 ! if(fftn2_distrib(i2)/=me_fft) cycle ! this ind is not to be treated by me_fft 2057 do i1=1,n1 2058 2059 ind1=ind1+1 2060 ! r2=ffti2_local(i2+1) - 1 2061 ! ind=n1*(nd2*i3+r2)+i1+1 !this is ind in the current proc 2062 2063 ! Get location of G vector (grid point) centered at 0 0 0 2064 l1=i1-(i1/id1)*n1-1 2065 l2=i2-(i2/id2)*n2-1 2066 l3=i3-(i3/id3)*n3-1 2067 2068 ! Get rotated G vector Gj for each symmetry element 2069 ! -- here we use the TRANSPOSE of symrel1; assuming symrel1 expresses 2070 ! the rotation in real space, the transpose is then appropriate 2071 ! for G space symmetrization (p. 1172d,e of notes, 2 June 1995). 2072 j1=symrel1(1,1)*l1+symrel1(2,1)*l2+symrel1(3,1)*l3 2073 j2=symrel1(1,2)*l1+symrel1(2,2)*l2+symrel1(3,2)*l3 2074 j3=symrel1(1,3)*l1+symrel1(2,3)*l2+symrel1(3,3)*l3 2075 2076 ! Map into [0,n-1] and then add 1 for array index in [1,n] 2077 k1=1+mod(n1+mod(j1,n1),n1) 2078 k2=1+mod(n2+mod(j2,n2),n2) 2079 k3=1+mod(n3+mod(j3,n3),n3) 2080 2081 ! Get linear index of rotated point Gj 2082 ind2=k1+n1*((k2-1)+n2*(k3-1)) 2083 ! r2=ffti2_local(j2+1) - 1 2084 ! ind=n1*(nd2*j3+r2)+j1+1 !this is ind may be in another proc!! 2085 2086 phnon1(1) = one 2087 phnon1(2) = zero 2088 if (t_tnon_nonzero) then 2089 ! compute exp(-2*Pi*I*G dot tau) using original G 2090 ! NB: this phase is same as that in irrzg and phnons1, and corresponds to complex conjugate of phase from G to Gj; 2091 ! we use it immediately below, to go _to_ workg(ind1) 2092 ! TODO : replace this with complex powers of exp(2pi tnon(1)) etc... 2093 arg=two_pi*(dble(l1)*tnon(1)+dble(l2)*tnon(2)+dble(l3)*tnon(3)) 2094 phnon1(1) = cos(arg) 2095 phnon1(2) =-sin(arg) 2096 end if 2097 2098 ! rho(Strans*G)=exp(2*Pi*I*(G) dot tau_S) rho(G) 2099 workg_eq (1, ind1) = phnon1(1) * workg(1, ind2) & 2100 & - phnon1(2) * workg(2, ind2) 2101 workg_eq (2, ind1) = phnon1(1) * workg(2, ind2) & 2102 & + phnon1(2) * workg(1, ind2) 2103 2104 end do 2105 end do 2106 end do 2107 2108 ! accumulate rhog1_eq 2109 if (ispden == 1) rhog1_eq = workg_eq 2110 2111 ! FFT back to real space to get rhor1_eq 2112 ! Pull out full or spin up density, now symmetrized 2113 call fourdp(cplex,workg_eq,rhor1_eq(:,ispden),1,mpi_enreg,nfft,ngfft,mpi_enreg%paral_kgb,0) 2114 2115 end do !nspden 2116 2117 ABI_DEALLOCATE(workg) 2118 ABI_DEALLOCATE(workg_eq) 2119 2120 end subroutine rotate_rho
m_spacepar/setsym [ Functions ]
[ Top ] [ m_spacepar ] [ Functions ]
NAME
setsym
FUNCTION
Set up irreducible zone in G space by direct calculation. Do not call this routine if nsym=1 (only identity symmetry). Only indsym and symrec get returned if iscf=0. symrec needed to symmetrize coordinate gradients in sygrad. (symrec is redundant and could be removed later in favor of symrel)
INPUTS
iscf=(<= 0 =>non-SCF), >0 => SCF natom=number of atoms in unit cell nfft=(effective) number of FFT grid points (for this processor) ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft nspden=number of spin-density components nsppol=1 for unpolarized, 2 for spin-polarized nsym=number of symmetries in space group (at least 1) symafm(nsym)=(anti)ferromagnetic part of symmetry operations symrel(3,3,nsym)=symmetry operations in terms of real space primitive translations tnons(3,nsym)=nonsymmorphic translations of space group in terms of real space primitive translations (may be 0) typat(natom)=atom type (integer) for each atom xred(3,natom)=atomic coordinates in terms of real space primitive translations
OUTPUT
indsym(4,nsym,natom)=indirect indexing of atom labels--see subroutine symatm for definition (if nsym>1) irrzon(nfft,2,(nspden/nsppol)-3*(nspden/4))=irreducible zone data phnons(2,nfft,(nspden/nsppol)-3*(nspden/4))=nonsymmorphic translation phases symrec(3,3,nsym)=symmetry operations in terms of reciprocal space primitive translations (if nsym>1)
NOTES
nsppol and nspden are needed in case of (anti)ferromagnetic symmetry operations
PARENTS
dfpt_looppert,gstate,m_dvdb,nonlinear,respfn,scfcv
CHILDREN
chkgrp,irrzg,mati3inv,symatm,symdet,timab
SOURCE
2168 subroutine setsym(indsym,irrzon,iscf,natom,nfft,ngfft,nspden,nsppol,nsym,phnons,& 2169 & symafm,symrec,symrel,tnons,typat,xred) 2170 2171 2172 !This section has been created automatically by the script Abilint (TD). 2173 !Do not modify the following lines by hand. 2174 #undef ABI_FUNC 2175 #define ABI_FUNC 'setsym' 2176 !End of the abilint section 2177 2178 implicit none 2179 2180 !Arguments ------------------------------------ 2181 !scalars 2182 integer,intent(in) :: iscf,natom,nfft,nspden,nsppol,nsym 2183 !arrays 2184 integer,intent(in) :: ngfft(18),symafm(nsym),symrel(3,3,nsym),typat(natom) 2185 integer,intent(out) :: indsym(4,nsym,natom) 2186 integer,intent(inout) :: irrzon(nfft,2,(nspden/nsppol)-3*(nspden/4)) !vz_i 2187 integer,intent(out) :: symrec(3,3,nsym) 2188 real(dp),intent(in) :: tnons(3,nsym),xred(3,natom) 2189 real(dp),intent(out) :: phnons(2,nfft,(nspden/nsppol)-3*(nspden/4)) 2190 2191 !Local variables------------------------------- 2192 !scalars 2193 integer :: isym,ierr 2194 real(dp) :: tolsym8 2195 !arrays 2196 integer,allocatable :: determinant(:) 2197 real(dp) :: tsec(2) 2198 2199 ! ************************************************************************* 2200 2201 call timab(6,1,tsec) 2202 2203 !Check that symmetries have unity determinant 2204 ABI_ALLOCATE(determinant,(nsym)) 2205 call symdet(determinant,nsym,symrel) 2206 ABI_DEALLOCATE(determinant) 2207 2208 2209 !Get the symmetry matrices in terms of reciprocal basis 2210 do isym=1,nsym 2211 call mati3inv(symrel(:,:,isym),symrec(:,:,isym)) 2212 end do 2213 2214 !Check for group closure 2215 call chkgrp(nsym,symafm,symrel,ierr) 2216 ABI_CHECK(ierr==0,"Error in group closure") 2217 2218 call chkgrp(nsym,symafm,symrec,ierr) 2219 ABI_CHECK(ierr==0,"Error in group closure") 2220 2221 !Obtain a list of rotated atom labels: 2222 tolsym8=tol8 2223 call symatm(indsym,natom,nsym,symrec,tnons,tolsym8,typat,xred) 2224 2225 !If non-SCF calculation, or nsym==1, do not need IBZ data 2226 if ( (iscf>0 .or. iscf==-3) .and. nsym>1 ) then 2227 ! Locate irreducible zone in reciprocal space for symmetrization: 2228 call irrzg(irrzon,nspden,nsppol,nsym,ngfft(1),ngfft(2),ngfft(3),phnons,symafm,symrel,tnons) 2229 end if 2230 2231 call timab(6,2,tsec) 2232 2233 end subroutine setsym
m_spacepar/symrhg [ Functions ]
[ Top ] [ m_spacepar ] [ Functions ]
NAME
symrhg
FUNCTION
From rho(r), generate rho(G), symmetrize it, and come back to the real space for a symmetrized rho(r).
INPUTS
cplex=1 if rhor is real, 2 if rhor is complex gprimd(3,3)=dimensional reciprocal space primitive translations irrzon(nfft,2,(nspden/nsppol)-3*(nspden/4))=irreducible zone data mpi_enreg=information about MPI parallelization nfft=(effective) number of FFT grid points (for this processor) ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft nspden=number of spin-density components nsppol=1 for unpolarized, 2 for spin-polarized nsym=number of symmetry elements. phnons(2,nfft,(nspden/nsppol)-3*(nspden/4))=nonsymmorphic translation phases rprimd(3,3)=dimensional real space primitive translations symafm(nsym)=(anti)ferromagnetic part of symmetry operations symrel(3,3,nsym)=symmetry matrices in real space (integers)
OUTPUT
rhog(2,nfft)=symmetrized rho(G) (total) electron density in G space
SIDE EFFECTS
Input/Output rhor(cplex*nfft,nspden)=array for electron density in electrons/bohr**3. Input, but also output, if symmetrization is applied. Also output if nspden > 1 (change spin components)
NOTES
When using spin-polarization (nspden==2), put total density in first half of rhor array and spin up in second half If (nspden=2 and nsppol=2) the density is transformed as (up,down) => (up+down,up) If (nspden=2 and nsppol=1) anti-ferromagnetic symmetry operations must be used, such as to transform (2*up) => (up+down,up) In spin-polarized, and if there is no symmetry to be applied on the system, only the total density is generated in G space
PARENTS
dfpt_mkrho,dfpt_nstpaw,dfpt_rhofermi,dfpt_vtorho,m_dvdb,mkrho suscep_stat,vtorho,vtorhorec,vtorhotf,wfd_mkrho
CHILDREN
fourdp,matr3inv,ptabs_fourdp,symredcart,timab,xmpi_sum
SOURCE
1116 subroutine symrhg(cplex,gprimd,irrzon,mpi_enreg,nfft,nfftot,ngfft,nspden,nsppol,nsym,paral_kgb,& 1117 & phnons,rhog,rhor,rprimd,symafm,symrel) 1118 1119 1120 !This section has been created automatically by the script Abilint (TD). 1121 !Do not modify the following lines by hand. 1122 #undef ABI_FUNC 1123 #define ABI_FUNC 'symrhg' 1124 !End of the abilint section 1125 1126 implicit none 1127 1128 !Arguments ------------------------------------ 1129 !scalars 1130 integer,intent(in) :: cplex,nfft,nfftot,nspden,nsppol,nsym,paral_kgb 1131 type(MPI_type),intent(in) :: mpi_enreg 1132 !arrays 1133 integer,intent(in) :: irrzon(nfftot**(1-1/nsym),2,(nspden/nsppol)-3*(nspden/4)),ngfft(18) 1134 integer,intent(in) :: symafm(nsym),symrel(3,3,nsym) 1135 real(dp),intent(in) :: gprimd(3,3),phnons(2,nfftot**(1-1/nsym),(nspden/nsppol)-3*(nspden/4)),rprimd(3,3) 1136 real(dp),intent(inout) :: rhor(cplex*nfft,nspden) 1137 real(dp),intent(out) :: rhog(2,nfft) 1138 1139 !Local variables------------------------------- 1140 !scalars 1141 integer :: ier,imagn,ind,ind2,indsy,ispden,isym,iup,izone,izone_max,j,j1,j2,j3,jsym 1142 integer :: k1,k2,k3,l1,l2,l3,me_fft 1143 integer :: n1,n2,n3,nd2,nproc_fft,nspden_eff,nsym_used,numpt,nup 1144 integer :: r2,rep,spaceComm 1145 logical,parameter :: afm_noncoll=.true. ! TRUE if antiferro symmetries are used in non-collinear magnetism 1146 real(dp) :: magxsu1,magxsu2,magysu1,magysu2,magzsu1,magzsu2,mxi,mxr,myi,myr,mzi,mzr,phi,phr,rhosu1,rhosu2 1147 !character(len=500) :: message 1148 !arrays 1149 integer,allocatable :: isymg(:) 1150 integer, ABI_CONTIGUOUS pointer :: fftn2_distrib(:),ffti2_local(:) 1151 integer, ABI_CONTIGUOUS pointer :: fftn3_distrib(:),ffti3_local(:) 1152 real(dp) :: tsec(2) 1153 real(dp),allocatable :: magngx(:,:),magngy(:,:),magngz(:,:) 1154 real(dp),allocatable :: rhosu1_arr(:),rhosu2_arr(:),work(:) 1155 real(dp),allocatable :: symafm_used(:),symrec_cart(:,:,:),symrel_cart(:,:,:) 1156 1157 !************************************************************************* 1158 ! 1159 !Note the timing channel 17 excludes the different Fourier transforms 1160 1161 ABI_ALLOCATE(work,(cplex*nfft)) 1162 1163 !Special treatment for spin-polarized case 1164 if(nspden==2 .and. nsppol==2) then 1165 ! When nspden=2 and nsppol=2, put total density in first half 1166 ! of rhor array and spin up in second half (up,down) => (up+down,up) 1167 call timab(17,1,tsec) 1168 work(:)=rhor(:,1) ! up => work 1169 rhor(:,1)=rhor(:,1)+rhor(:,2) ! up+down 1170 rhor(:,2)=work(:) ! work => up 1171 call timab(17,2,tsec) 1172 end if 1173 1174 !Special treatment for antiferromagnetism case 1175 if(nspden==2 .and. nsppol==1) then 1176 call timab(17,1,tsec) 1177 ! When nspden=2 and nsppol=1, (2*up) => (2*up,up) 1178 ! Indeed, what was delivered to the present routine is a "total" density, 1179 ! obtained from occupation numbers varying between 0 and 2, 1180 ! but for spin up only potential. 1181 rhor(:,2)=half*rhor(:,1) 1182 call timab(17,2,tsec) 1183 end if 1184 1185 !Special treatment for non-collinear magnetism case 1186 if(nspden==4) then 1187 call timab(17,1,tsec) 1188 !FR the half factors missing are recovered in dfpt_mkvxc_noncoll and dfpt_accrho 1189 rhor(:,1)=rhor(:,1)+rhor(:,4) !nup+ndown 1190 rhor(:,2)=rhor(:,2)-rhor(:,1) !mx (n+mx-n) 1191 rhor(:,3)=rhor(:,3)-rhor(:,1) !my (n+my-n) 1192 rhor(:,4)=rhor(:,1)-two*rhor(:,4) !mz=n-2ndown 1193 call timab(17,2,tsec) 1194 end if 1195 1196 1197 if(nsym==1)then 1198 1199 if(nspden==2 .and. nsppol==1) then ! There must be at least one anti-ferromagnetic operation 1200 MSG_BUG('In the antiferromagnetic case, nsym cannot be 1') 1201 end if 1202 1203 ! If not using symmetry, still want total density in G space rho(G). 1204 ! Fourier transform (incl normalization) to get rho(G) 1205 work(:)=rhor(:,1) 1206 call fourdp(cplex,rhog,work,-1,mpi_enreg,nfft,ngfft,paral_kgb,0) 1207 else 1208 1209 ! Treat either full density, spin-up density or magnetization 1210 ! Note the decrease of ispden to the value 1, in order to finish 1211 ! with rhog of the total density (and not the spin-up density or magnetization) 1212 nspden_eff=nspden;if (nspden==4) nspden_eff=1 1213 do ispden=nspden_eff,1,-1 1214 1215 ! Prepare the density to be symmetrized, in the reciprocal space 1216 if(nspden==1 .or. nsppol==2 .or. (nspden==4.and.(.not.afm_noncoll)))then 1217 imagn=1 1218 nsym_used=0 1219 do isym=1,nsym 1220 if(symafm(isym)==1)nsym_used=nsym_used+1 1221 ! DEBUG 1222 ! write(std_out,*)' symrhg : isym,symafm(isym)',isym,symafm(isym) 1223 ! ENDDEBUG 1224 end do 1225 else if(nspden==2 .and. nsppol==1)then ! antiferromagnetic case 1226 imagn=ispden 1227 nsym_used=nsym/ispden 1228 else if (nspden==4) then 1229 imagn=1 1230 nsym_used=nsym/ispden 1231 end if 1232 1233 ! write(std_out,*)' symrhg : nsym_used=',nsym_used 1234 1235 ! rhor -fft-> rhog (rhog is used as work space) 1236 ! Note : it should be possible to reuse rhog in the antiferromagnetic case this would avoid one FFT 1237 work(:)=rhor(:,ispden) 1238 call fourdp(cplex,rhog,work,-1,mpi_enreg,nfft,ngfft,paral_kgb,0) 1239 if (nspden==4) then 1240 ABI_ALLOCATE(magngx,(2,nfft)) 1241 ABI_ALLOCATE(magngy,(2,nfft)) 1242 ABI_ALLOCATE(magngz,(2,nfft)) 1243 work(:)=rhor(:,2) 1244 call fourdp(cplex,magngx,work,-1,mpi_enreg,nfft,ngfft,paral_kgb,0) 1245 work(:)=rhor(:,3) 1246 call fourdp(cplex,magngy,work,-1,mpi_enreg,nfft,ngfft,paral_kgb,0) 1247 work(:)=rhor(:,4) 1248 call fourdp(cplex,magngz,work,-1,mpi_enreg,nfft,ngfft,paral_kgb,0) 1249 end if 1250 1251 ! Begins the timing here only , to exclude FFTs 1252 call timab(17,1,tsec) 1253 1254 n1=ngfft(1);n2=ngfft(2);n3=ngfft(3);nproc_fft=ngfft(10);me_fft=ngfft(11);nd2=n2/nproc_fft 1255 1256 ! Get the distrib associated with this fft_grid 1257 call ptabs_fourdp(mpi_enreg,n2,n3,fftn2_distrib,ffti2_local,fftn3_distrib,ffti3_local) 1258 1259 ! The following is only valid for total, up or dn density 1260 ! ------------------------------------------------------- 1261 1262 ! Get maxvalue of izone 1263 izone_max=count(irrzon(:,2,imagn)>0) 1264 ABI_ALLOCATE(rhosu1_arr,(izone_max)) 1265 ABI_ALLOCATE(rhosu2_arr,(izone_max)) 1266 1267 numpt=0 1268 do izone=1,nfftot 1269 1270 ! Get repetition number 1271 rep=irrzon(izone,2,imagn) 1272 if(rep==0)exit 1273 1274 ! Compute number of unique points in this symm class: 1275 nup=nsym_used/rep 1276 1277 ! Accumulate charge over equivalent points 1278 rhosu1=zero 1279 rhosu2=zero 1280 do iup=1,nup 1281 ind=irrzon(iup+numpt,1,imagn) 1282 j=ind-1;j1=modulo(j,n1);j2=modulo(j/n1,n2);j3=j/(n1*n2); 1283 if(fftn2_distrib(j2+1)==me_fft) then ! this ind is to be treated by me_fft 1284 r2=ffti2_local(j2+1) - 1 1285 ind=n1*(nd2*j3+r2)+j1+1 !this is ind in the current proc 1286 rhosu1=rhosu1+rhog(1,ind)*phnons(1,iup+numpt,imagn)& 1287 & -rhog(2,ind)*phnons(2,iup+numpt,imagn) 1288 rhosu2=rhosu2+rhog(2,ind)*phnons(1,iup+numpt,imagn)& 1289 & +rhog(1,ind)*phnons(2,iup+numpt,imagn) 1290 end if 1291 1292 end do 1293 rhosu1=rhosu1/dble(nup) 1294 rhosu2=rhosu2/dble(nup) 1295 rhosu1_arr(izone)=rhosu1 1296 rhosu2_arr(izone)=rhosu2 1297 ! Keep index of how many points have been considered: 1298 numpt=numpt+nup 1299 1300 end do ! End loop over izone 1301 1302 ! Reduction in case of FFT parallelization 1303 if(mpi_enreg%nproc_fft>1)then 1304 spaceComm=mpi_enreg%comm_fft 1305 call xmpi_sum(rhosu1_arr,spaceComm,ier) 1306 call xmpi_sum(rhosu2_arr,spaceComm,ier) 1307 end if 1308 1309 ! Now symmetrize the density 1310 numpt=0 1311 do izone=1,nfftot 1312 1313 ! Get repetition number 1314 rep=irrzon(izone,2,imagn) 1315 if(rep==0)exit 1316 1317 ! Compute number of unique points in this symm class: 1318 nup=nsym_used/rep 1319 1320 ! Define symmetrized rho(G) at equivalent points: 1321 do iup=1,nup 1322 ind=irrzon(iup+numpt,1,imagn) 1323 ! decompose ind-1=n1(n2 j3+ j2)+j1 1324 j=ind-1;j1=modulo(j,n1);j2=modulo(j/n1,n2);j3=j/(n1*n2); 1325 if(fftn2_distrib(j2+1)==me_fft) then ! this ind is to be treated by me_fft 1326 r2=ffti2_local(j2+1) - 1 1327 ! ind in the proc ind-1=n1(nd2 j3+ r2)+j1 1328 ind=n1*(nd2*j3+r2)+j1+1 !this is ind in the current proc 1329 rhog(1,ind)=rhosu1_arr(izone)*phnons(1,iup+numpt,imagn)& 1330 & +rhosu2_arr(izone)*phnons(2,iup+numpt,imagn) 1331 rhog(2,ind)=rhosu2_arr(izone)*phnons(1,iup+numpt,imagn)& 1332 & -rhosu1_arr(izone)*phnons(2,iup+numpt,imagn) 1333 end if 1334 end do 1335 1336 ! Keep index of how many points have been considered: 1337 numpt=numpt+nup 1338 1339 end do ! End loop over izone 1340 1341 ABI_DEALLOCATE(rhosu1_arr) 1342 ABI_DEALLOCATE(rhosu2_arr) 1343 1344 ! The following is only valid for magnetization 1345 ! --------------------------------------------- 1346 if (nspden==4) then 1347 1348 ! Transfer symmetries in cartesian coordinates 1349 ! Compute symmetries in reciprocal space in cartesian coordinates 1350 ABI_ALLOCATE(symrec_cart,(3,3,nsym_used)) 1351 ABI_ALLOCATE(symrel_cart,(3,3,nsym_used)) 1352 ABI_ALLOCATE(symafm_used,(nsym_used)) 1353 jsym=0 1354 do isym=1,nsym 1355 if (symafm(isym)/=1.and.(.not.afm_noncoll)) cycle 1356 jsym=jsym+1 1357 symafm_used(jsym)=dble(symafm(isym)) 1358 call symredcart(rprimd,gprimd,symrel_cart(:,:,jsym),symrel(:,:,isym)) 1359 call matr3inv(symrel_cart(:,:,jsym),symrec_cart(:,:,jsym)) 1360 end do 1361 1362 numpt=count(irrzon(:,1,imagn)>0) 1363 ABI_ALLOCATE(isymg,(numpt)) 1364 isymg=0 1365 ABI_ALLOCATE(rhosu1_arr,(3*izone_max)) 1366 ABI_ALLOCATE(rhosu2_arr,(3*izone_max)) 1367 1368 ! Accumulate magnetization over equivalent points 1369 ! Use all symmetries (not only those linking different g points) 1370 ! Use Inverse[Transpose[symrel]]=symrec 1371 numpt=0 1372 do izone=1,izone_max 1373 magxsu1=zero;magxsu2=zero 1374 magysu1=zero;magysu2=zero 1375 magzsu1=zero;magzsu2=zero 1376 ind=irrzon(1+numpt,1,1) 1377 rep=irrzon(izone,2,1) 1378 nup=nsym_used/rep 1379 j=ind-1;l1=modulo(j,n1);l2=modulo(j/n1,n2);l3=j/(n1*n2) 1380 jsym=0 1381 do isym=1,nsym 1382 if (symafm(isym)/=1.and.(.not.afm_noncoll)) cycle 1383 jsym=jsym+1 1384 j1=symrel(1,1,isym)*l1+symrel(2,1,isym)*l2+symrel(3,1,isym)*l3 1385 j2=symrel(1,2,isym)*l1+symrel(2,2,isym)*l2+symrel(3,2,isym)*l3 1386 j3=symrel(1,3,isym)*l1+symrel(2,3,isym)*l2+symrel(3,3,isym)*l3 1387 k1=map_symrhg(j1,n1);k2=map_symrhg(j2,n2);k3=map_symrhg(j3,n3) 1388 indsy=1+k1+n1*(k2+n2*k3) 1389 ind2=-1;iup=numpt 1390 do while (ind2/=indsy.and.iup<numpt+nup) 1391 iup=iup+1;ind2=irrzon(iup,1,1) 1392 end do 1393 if (ind2/=indsy) then 1394 MSG_ERROR("ind2/=indsy in symrhg !") 1395 end if 1396 if (isymg(iup)==0) isymg(iup)=jsym 1397 if(fftn2_distrib(modulo((indsy-1)/n1,n2) + 1) == me_fft ) then ! this is indsy is to be treated by me_fft 1398 indsy=n1*(nd2*k3+ ffti2_local(k2+1) -1)+k1+1 ! this is indsy in the current proc 1399 phr=phnons(1,iup,imagn);if (rep==1) phr=phr*symafm_used(jsym) !if rep==2, symafm is already included in phnons 1400 phi=phnons(2,iup,imagn);if (rep==1) phi=phi*symafm_used(jsym) !(see irrzg.F90) 1401 mxr=symrel_cart(1,1,jsym)*magngx(1,indsy)+symrel_cart(1,2,jsym)*magngy(1,indsy)+symrel_cart(1,3,jsym)*magngz(1,indsy) 1402 mxi=symrel_cart(1,1,jsym)*magngx(2,indsy)+symrel_cart(1,2,jsym)*magngy(2,indsy)+symrel_cart(1,3,jsym)*magngz(2,indsy) 1403 myr=symrel_cart(2,1,jsym)*magngx(1,indsy)+symrel_cart(2,2,jsym)*magngy(1,indsy)+symrel_cart(2,3,jsym)*magngz(1,indsy) 1404 myi=symrel_cart(2,1,jsym)*magngx(2,indsy)+symrel_cart(2,2,jsym)*magngy(2,indsy)+symrel_cart(2,3,jsym)*magngz(2,indsy) 1405 mzr=symrel_cart(3,1,jsym)*magngx(1,indsy)+symrel_cart(3,2,jsym)*magngy(1,indsy)+symrel_cart(3,3,jsym)*magngz(1,indsy) 1406 mzi=symrel_cart(3,1,jsym)*magngx(2,indsy)+symrel_cart(3,2,jsym)*magngy(2,indsy)+symrel_cart(3,3,jsym)*magngz(2,indsy) 1407 magxsu1=magxsu1+mxr*phr-mxi*phi;magxsu2=magxsu2+mxi*phr+mxr*phi 1408 magysu1=magysu1+myr*phr-myi*phi;magysu2=magysu2+myi*phr+myr*phi 1409 magzsu1=magzsu1+mzr*phr-mzi*phi;magzsu2=magzsu2+mzi*phr+mzr*phi 1410 end if 1411 end do 1412 rhosu1_arr(3*izone-2)=magxsu1/dble(nsym_used) 1413 rhosu1_arr(3*izone-1)=magysu1/dble(nsym_used) 1414 rhosu1_arr(3*izone )=magzsu1/dble(nsym_used) 1415 rhosu2_arr(3*izone-2)=magxsu2/dble(nsym_used) 1416 rhosu2_arr(3*izone-1)=magysu2/dble(nsym_used) 1417 rhosu2_arr(3*izone )=magzsu2/dble(nsym_used) 1418 numpt=numpt+nup 1419 end do 1420 1421 ! Reduction in case of FFT parallelization 1422 if(mpi_enreg%nproc_fft>1)then 1423 spaceComm=mpi_enreg%comm_fft 1424 call xmpi_sum(rhosu1_arr,spaceComm,ier) 1425 call xmpi_sum(rhosu2_arr,spaceComm,ier) 1426 end if 1427 1428 ! Now symmetrize the magnetization at equivalent points 1429 ! Use Transpose[symrel] 1430 numpt=0 1431 do izone=1,izone_max 1432 rep=irrzon(izone,2,imagn) 1433 nup=nsym_used/rep 1434 do iup=1,nup 1435 ind=irrzon(iup+numpt,1,imagn) 1436 j=ind-1;j1=modulo(j,n1);j2=modulo(j/n1,n2);j3=j/(n1*n2) 1437 if(fftn2_distrib(j2+1)==me_fft) then ! this ind is to be treated by me_fft 1438 r2=ffti2_local(j2+1) - 1 1439 ind=n1*(nd2*j3+r2)+j1+1 ! this is ind in the current proc 1440 jsym=isymg(iup+numpt) 1441 if (jsym==0) then 1442 MSG_ERROR("jsym=0 in symrhg !") 1443 end if 1444 magxsu1=rhosu1_arr(3*izone-2);magxsu2=rhosu2_arr(3*izone-2) 1445 magysu1=rhosu1_arr(3*izone-1);magysu2=rhosu2_arr(3*izone-1) 1446 magzsu1=rhosu1_arr(3*izone );magzsu2=rhosu2_arr(3*izone ) 1447 phr=phnons(1,iup,imagn);if (rep==1) phr=phr*symafm_used(jsym) !if rep==2, symafm is already included in phnons 1448 phi=phnons(2,iup,imagn);if (rep==1) phi=phi*symafm_used(jsym) !(see irrzg.F90) 1449 mxr=symrec_cart(1,1,jsym)*magxsu1+symrec_cart(2,1,jsym)*magysu1+symrec_cart(3,1,jsym)*magzsu1 1450 mxi=symrec_cart(1,1,jsym)*magxsu2+symrec_cart(2,1,jsym)*magysu2+symrec_cart(3,1,jsym)*magzsu2 1451 myr=symrec_cart(1,2,jsym)*magxsu1+symrec_cart(2,2,jsym)*magysu1+symrec_cart(3,2,jsym)*magzsu1 1452 myi=symrec_cart(1,2,jsym)*magxsu2+symrec_cart(2,2,jsym)*magysu2+symrec_cart(3,2,jsym)*magzsu2 1453 mzr=symrec_cart(1,3,jsym)*magxsu1+symrec_cart(2,3,jsym)*magysu1+symrec_cart(3,3,jsym)*magzsu1 1454 mzi=symrec_cart(1,3,jsym)*magxsu2+symrec_cart(2,3,jsym)*magysu2+symrec_cart(3,3,jsym)*magzsu2 1455 magngx(1,ind)=mxr*phr+mxi*phi 1456 magngx(2,ind)=mxi*phr-mxr*phi 1457 magngy(1,ind)=myr*phr+myi*phi 1458 magngy(2,ind)=myi*phr-myr*phi 1459 magngz(1,ind)=mzr*phr+mzi*phi 1460 magngz(2,ind)=mzi*phr-mzr*phi 1461 end if 1462 end do 1463 numpt=numpt+nup 1464 end do 1465 ABI_DEALLOCATE(isymg) 1466 ABI_DEALLOCATE(rhosu1_arr) 1467 ABI_DEALLOCATE(rhosu2_arr) 1468 ABI_DEALLOCATE(symrec_cart) 1469 ABI_DEALLOCATE(symrel_cart) 1470 ABI_DEALLOCATE(symafm_used) 1471 1472 end if ! nspden==4 1473 1474 call timab(17,2,tsec) 1475 1476 ! Pull out full or spin up density, now symmetrized 1477 call fourdp(cplex,rhog,work,1,mpi_enreg,nfft,ngfft,paral_kgb,0) 1478 rhor(:,ispden)=work(:) 1479 if (nspden==4) then 1480 call fourdp(cplex,magngx,work,1,mpi_enreg,nfft,ngfft,paral_kgb,0) 1481 rhor(:,2)=work(:) 1482 call fourdp(cplex,magngy,work,1,mpi_enreg,nfft,ngfft,paral_kgb,0) 1483 rhor(:,3)=work(:) 1484 call fourdp(cplex,magngz,work,1,mpi_enreg,nfft,ngfft,paral_kgb,0) 1485 rhor(:,4)=work(:) 1486 ABI_DEALLOCATE(magngx) 1487 ABI_DEALLOCATE(magngy) 1488 ABI_DEALLOCATE(magngz) 1489 end if 1490 1491 end do ! ispden 1492 1493 end if ! End on the condition nsym==1 1494 1495 ABI_DEALLOCATE(work) 1496 1497 contains 1498 1499 function map_symrhg(j1,n1) 1500 1501 1502 !This section has been created automatically by the script Abilint (TD). 1503 !Do not modify the following lines by hand. 1504 #undef ABI_FUNC 1505 #define ABI_FUNC 'map_symrhg' 1506 !End of the abilint section 1507 1508 integer :: map_symrhg 1509 integer,intent(in) :: j1,n1 1510 map_symrhg=mod(n1+mod(j1,n1),n1) 1511 end function map_symrhg 1512 1513 end subroutine symrhg