TABLE OF CONTENTS
ABINIT/wfconv [ Functions ]
NAME
wfconv
FUNCTION
This subroutine treats the wavefunctions for one k point, and converts them to other parameters.
COPYRIGHT
Copyright (C) 2000-2018 ABINIT group (XG,TD) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt . For the initials of contributors, see ~abinit/doc/developers/contributors.txt .
INPUTS
ceksp2=if 1, center the output sphere of pw on Gamma; if 0, on each k-point (usual). cg1(2,mcg1)=wavefunction array debug= if 1, print some messages ; otherwise, 0. ecut1=kinetic energy cutoffs for basis sphere 1 (hartree) ecut2=kinetic energy cutoff beyond which the coefficients of wf2 vanish (Ha) ecut2_eff=kinetic energy cut-off for basis sphere 2 (hartree) eig_k1(mband1*(2*mband1)**formeig)=eigenvalues exchn2n3d=if 1, n2 and n3 are exchanged formeig option (format of the eigenvalues and eigenvector) : 0 => ground-state format (initialisation of eigenvectors with random numbers, vector of eigenvalues) 1 => respfn format (initialisation of eigenvectors with 0 s, hermitian matrix of eigenvalues) gmet1(3,3)=reciprocal space metric (bohr^-2) for input wf gmet2(3,3)=reciprocal space metric (bohr^-2) for output wf icg1=shift to be given to the location of the data in the array cg1 icg2=shift to be given to the location of the data in the array cg2 ikpt1=number of the k point actually treated (input wf numbering) ikpt10=number of the k point previously treated (input wf numbering) ikpt2=number of the k point actually treated (output numbering) indkk(nkpt2*sppoldbl,6)=describe k point number of kptns1 that allows to generate wavefunctions closest to given kpt2 (and possibly isppol2=2) indkk(:,1)=k point number of kpt1 indkk(:,2)=symmetry operation to be applied to kpt1, to give kpt1a (if 0, means no symmetry operation, equivalent to identity ) indkk(:,3:5)=shift in reciprocal space to be given to kpt1a, to give kpt1b, that is the closest to kpt2. indkk(:,6)=1 if time-reversal was used to generate kpt1a from kpt1, 0 otherwise inplace= if 0, cg1 and cg2 are different in the calling routine, if 1, cg1 and cg2 are identical (they have the same memory location) This is also true for the pairs (eig_k1,eig_k2) and (occ_k1,occ_k2) isppol2=spin variable for output wavefunctions istwfk1(nkpt1)=input parameter that describes the storage of wfs in set1 istwfk2(nkpt2)=input parameter that describes the storage of wfs in set2 kg1(3,mpw1)=dimensionless coords of G vecs in basis sphere at k point (input wf) kg2(3,mpw2)=dimensionless coords of G vecs in basis sphere at k point (output wf) kptns1(3,nkpt1)=k point set for input wavefunctions kptns2(3,nkpt2)=k point set for output wavefunctions mband1=dimension of eig_k1 and occ_k1 arrays mband2=dimension of eig_k2 and occ_k2 arrays mcg1=dimension of cg1 array (at least npw1*nspinor1*nbd1) mcg2=dimension of cg2 array (at least npw2*nspinor2*nbd2) mpi_enreg1=informations about MPI parallelization for set 1 mpi_enreg2=informations about MPI parallelization for set 2 mpw1=dimension of kg1, can be set to 0 if not needed mpw2=dimension of kg2, can be set to 0 if not needed nbd1=number of bands contained in cg1,eig_k1,occ_k1 at this k-point - spin (at input) nbd2=number of bands contained in cg2,eig_k2,occ_k2 at this k-point - spin (at output) ngfft1(18)=all needed information about 3D FFT, for input wavefunctions ngfft2(18)=all needed information about 3D FFT, for output wavefunctions see ~abinit/doc/variables/vargs.htm#ngfft nkpt1=number of k points for input wavefunctions nkpt2=number of k points for output wavefunctions npw1=number of planewaves for input wavefunctions npw2=number of planewaves for output wavefunctions nspinor1=number of spinors for input wavefunctions nspinor2=number of spinors for output wavefunctions nsym=number of symmetry elements in space group occ_k1(mband1)=occupation numbers optorth=1 if the WFs are orthogonalized before leaving the routine randalg=1 if "good" (but non-portable) random numbers should be used, 0 for compatibility restart=if 2, conversion between wavefunctions if 1, direct restart is allowed (see hdr_check.f) rprimd2(3,3)=dimensional primitive translations for real space (bohr) needed only for the spinor rotation sppoldbl= if 1, no doubling of the number if spins thanks to antiferromagn if 2, deduce nsppol=2 from nsppol=1, using Shubnikov symmetries symrel(3,3,nsym)=symmetry operations in real space in terms of primitive translations tnons(3,nsym)=nonsymmorphic translations for symmetry operations
OUTPUT
cg2(2,mcg2)=wavefunction array eig_k2(mband2*(2*mband2)**formeig)=eigenvalues occ_k2(mband2)=occupation (completed with zeros)
SIDE EFFECTS
Input/Output:
ikpt10=at input, number of the k point previously treated (input wf numbering)
(if this is the first call for the present k point set, ikpt10 should be 0)
at output, number of the k point just treated (input wf numbering)
kg1, kg2, npw1 and npw2 should not be modified by kpgsph (TD).
NOTES
Note that this routine can make an in-place conversion (see the input variable "inplace"), if cg1 and cg2 are equal, as well as the pairs (icg1,icg2), (eig_k1,eig_k2),(occ_k1,occ_k2) and (mband1,mband2) It can also be used to fill or to initialize wavefunctions at one k point (filling with random numbers or 0''s, according to the value of formeig), if the input number of bands (nbd1) is 0. In the latter case, one should use the same values of input wavefunction parameters than for output wavefunction parameters, except nbd1. The input parameters are indexed with 1, the output parameters are indexed with 2. Some of the arguments are arrays dimensioned with nkpt1 or nkpt2. Note that for these, only the elements for ikpt1 or ikpt2 will be used. The number of input bands must already be minimal at the input. This means, when input and output nspinor are equal : nbd1<nbd2 When the two nspinor differ, one must have nbd1/nspinor1<nbd2/nspinor2
PARENTS
newkpt,wfsinp
CHILDREN
cg_envlop,getph,getspinrot,kpgsph,mati3inv,ph1d3d,pw_orthon,sphere sphereboundary,timab,wrtout,xmpi_sum
SOURCE
134 #if defined HAVE_CONFIG_H 135 #include "config.h" 136 #endif 137 138 #include "abi_common.h" 139 140 subroutine wfconv(ceksp2,cg1,cg2,debug,ecut1,ecut2,ecut2_eff,& 141 & eig_k1,eig_k2,exchn2n3d,formeig,gmet1,gmet2,icg1,icg2,& 142 & ikpt1,ikpt10,ikpt2,indkk,inplace,isppol2,istwfk1,istwfk2,& 143 & kg1,kg2,kptns1,kptns2,mband1,mband2,mcg1,mcg2,mpi_enreg1,mpi_enreg2,& 144 & mpw1,mpw2,nbd1,nbd2,ngfft1,ngfft2,nkpt1,nkpt2,npw1,npw2,nspinor1,nspinor2,& 145 & nsym,occ_k1,occ_k2,optorth,randalg,restart,rprimd2,sppoldbl,symrel,tnons) 146 147 use defs_basis 148 use defs_abitypes 149 use m_errors 150 use m_profiling_abi 151 use m_xmpi 152 153 use m_fftcore, only : kpgsph, sphere 154 use m_cgtools, only : cg_envlop 155 use m_kg, only : ph1d3d, getph 156 157 !This section has been created automatically by the script Abilint (TD). 158 !Do not modify the following lines by hand. 159 #undef ABI_FUNC 160 #define ABI_FUNC 'wfconv' 161 use interfaces_14_hidewrite 162 use interfaces_18_timing 163 use interfaces_32_util 164 use interfaces_41_geometry 165 use interfaces_52_fft_mpi_noabirule 166 use interfaces_66_wfs, except_this_one => wfconv 167 !End of the abilint section 168 169 implicit none 170 171 !Arguments ------------------------------------ 172 !scalars 173 integer,intent(in) :: ceksp2,debug,exchn2n3d,formeig,icg1,icg2,ikpt1 174 integer,intent(in) :: ikpt2,inplace,isppol2,mband1,mband2,mcg1,mcg2,mpw1,mpw2 175 integer,intent(in) :: nbd1,nbd2,nkpt1,nkpt2,nspinor1,nspinor2,nsym 176 integer,intent(in) :: optorth,randalg,restart,sppoldbl 177 integer,intent(inout) :: ikpt10,npw1,npw2 178 real(dp),intent(in) :: ecut1,ecut2,ecut2_eff 179 type(MPI_type),intent(inout) :: mpi_enreg1,mpi_enreg2 180 !arrays 181 integer,intent(in) :: indkk(nkpt2*sppoldbl,6),istwfk1(nkpt1),istwfk2(nkpt2) 182 integer,intent(in) :: ngfft1(18),ngfft2(18),symrel(3,3,nsym) 183 integer,intent(inout) :: kg1(3,mpw1),kg2(3,mpw2) 184 real(dp),intent(in) :: gmet1(3,3),gmet2(3,3),kptns1(3,nkpt1),kptns2(3,nkpt2) 185 real(dp),intent(in) :: rprimd2(3,3),tnons(3,nsym) 186 real(dp),intent(inout) :: cg1(2,mcg1),cg2(2,mcg2) 187 real(dp),intent(inout) :: eig_k1(mband1*(2*mband1)**formeig) 188 real(dp),intent(inout) :: eig_k2(mband2*(2*mband2)**formeig),occ_k1(mband1) 189 real(dp),intent(inout) :: occ_k2(mband2) 190 191 !Local variables ------------------------------ 192 !scalars 193 integer,parameter :: nkpt_max=50,tobox=1,tosph=-1 194 integer :: conv_tnons,convert,fftalg,fold1,fold2,foldim,foldre,i1,i2,iband 195 integer :: iband_first,iband_last,icgmod,ierr,index,ipw 196 integer :: ispinor,ispinor1,ispinor2,ispinor_first,ispinor_last 197 integer :: istwf10_k,istwf1_k,istwf2_k,isym,itimrev,jsign 198 integer :: mgfft1,mgfft2,n1,n2,n3,n4,n5,n6 199 integer :: nbremn,npwtot,nspinor_index,nspinor1_this_proc,nspinor2_this_proc 200 integer :: order,ortalgo,seed 201 real(dp) :: ai,ar,arg,bi,br,eig_tmp,spinrots,spinrotx,spinroty,spinrotz 202 character(len=500) :: message 203 integer, parameter :: int64 = selected_int_kind(18) 204 !arrays 205 integer :: atindx(1),identity(3,3),ngfft_now(18),no_shift(3),shiftg(3) 206 integer :: symm(3,3),symrel_conv(3,3) 207 integer,allocatable :: gbound1(:,:),gbound2(:,:) 208 real(dp) :: kpoint1(3),kpoint2_sph(3),phktnons(2,1),spinrot(4),tnons_conv(3),tsec(2) 209 real(dp),allocatable :: cfft(:,:,:,:),dum(:,:),phase1d(:,:),phase3d(:,:) 210 real(dp),allocatable :: wavef1(:,:),wavef2(:,:),wavefspinor(:,:) 211 212 ! ************************************************************************* 213 214 mgfft1=maxval(ngfft1(1:3)) 215 mgfft2=maxval(ngfft2(1:3)) 216 if(.false.)write(std_out,*)occ_k1 ! just to keep occ_k1 as an argument before resolving the issue of its transfer 217 218 if(nspinor1/=1 .and. nspinor1/=2)then 219 write(message,'(a,i0)')'The argument nspinor1 must be 1 or 2, while it is nspinor1 = ',nspinor1 220 MSG_BUG(message) 221 end if 222 223 if(nspinor2/=1 .and. nspinor2/=2)then 224 write(message,'(a,i0)')' The argument nspinor2 must be 1 or 2, while it is nspinor2=',nspinor2 225 MSG_BUG(message) 226 end if 227 228 if(nspinor1==2 .and. mod(nbd1,2)/=0)then 229 write(message,'(a,i0)')' When nspinor1 is 2, nbd1 must be even, while it is nbd1 = ',nbd1 230 MSG_BUG(message) 231 end if 232 233 if(nspinor2==2 .and. mod(nbd2,2)/=0)then 234 write(message,'(a,i0)')' When nspinor2 is 2, nbd2 must be even, while it is nbd2=',nbd2 235 MSG_BUG(message) 236 end if 237 238 if(nbd1/nspinor1>nbd2/nspinor2)then 239 write(message, '(3a,2i6,3a,2i6,a)' )& 240 & 'In wfconv, the nbd/nspinor ratio cannot decrease. However,',ch10,& 241 & 'the initial quantities are nbd1,nspinor1=',nbd1,nspinor1,', and',ch10,& 242 & 'the requested final quantities are nbd2,nspinor2=',nbd2,nspinor2,'.' 243 MSG_BUG(message) 244 end if 245 246 ngfft_now(1:3)=ngfft1(1:3) 247 ngfft_now(8:18)=ngfft1(8:18) 248 !This line is the reason why ngfft_now has to be introduced 249 ngfft_now(7)=101 250 ngfft_now(4:6)=ngfft_now(1:3) 251 n1=ngfft_now(1) ; n2=ngfft_now(2) ; n3=ngfft_now(3) 252 n4=ngfft_now(4) ; n5=ngfft_now(5) ; n6=ngfft_now(6) 253 fftalg=ngfft_now(7) 254 255 !Parallelization over spinors management 256 nspinor1_this_proc=max(1,nspinor1/mpi_enreg1%nproc_spinor) 257 nspinor2_this_proc=max(1,nspinor2/mpi_enreg2%nproc_spinor) 258 259 !In order to generate IN PLACE new wfs from old wfs, the loop 260 !over bands and spinors must be done in one direction or the other, 261 !depending on npw1 and npw2, nspinor1 and nspinor2. 262 !If nspinor1=1 and nspinor2=2 , note that one will generate 263 !from nbd1 states of npw1 coefficients, 264 !2*nbd1 states of 2*npw2 coefficients. nbd1 cancels in comparing 265 !these expressions, but nspinor2 appears squared. 266 !The same line of thought works for the case nspinor1=2 and nspinor2=1 267 order=1 268 iband_first=1 ; iband_last=nbd1 269 ispinor_first=1 ; ispinor_last=nspinor1 270 if(nspinor1==2 .and. nspinor2==1)then 271 order=2 ; iband_last=nbd1-1 ; ispinor_last=1 272 end if 273 !Here, reverse the order if needed 274 if( npw2*nspinor2**2 > npw1*nspinor1**2 )then 275 order=-order 276 iband_first=iband_last ; iband_last=1 277 ispinor_first=ispinor_last ; ispinor_last=1 278 end if 279 280 kpoint1(:)=kptns1(:,ikpt1) 281 istwf1_k=istwfk1(ikpt1) 282 283 kpoint2_sph(:)=0.0_dp 284 if(ceksp2==0)kpoint2_sph(:)=kptns2(:,ikpt2) 285 istwf2_k=istwfk2(ikpt2) 286 287 !DEBUG 288 !write(std_out,*)'ecut1,ecut2_eff=',ecut1,ecut2_eff 289 !write(std_out,*)'gmet1,gmet2=',gmet1,gmet2 290 !write(std_out,*)'kpoint1,kpoint2_sph=',kpoint1,kpoint2_sph 291 !write(std_out,*)'nspinor1,nspinor2',nspinor1,nspinor2 292 !write(std_out,*)'istwf1_k,istwf2_k=',istwf1_k,istwf2_k 293 !write(std_out,*)'nbd1,tol8=',nbd1,tol8 294 !ENDDEBUG 295 296 !Determine whether it will be needed to convert the existing 297 !wavefunctions, or simply to complete them. 298 299 convert=0 300 if(nbd1/=0)then 301 if(abs(ecut2_eff-ecut1)>tol8)convert=convert+1 302 if(sum(abs(gmet2(:,:)-gmet1(:,:)))>tol8)convert=convert+2 303 if(sum(abs(kpoint2_sph(:)-kpoint1(:)))>tol8)convert=convert+4 304 if(nspinor2/=nspinor1)convert=convert+8 305 if(istwf2_k/=istwf1_k)convert=convert+16 306 end if 307 308 !This is a supplementary check 309 if(restart==1 .and. convert/=0)then 310 MSG_BUG('Restart==1 and convert/=0 are exclusive') 311 end if 312 313 !Determine whether symmetries must be used 314 conv_tnons=0 315 no_shift(:)=0 316 identity(:,:)=0 317 identity(1,1)=1 ; identity(2,2)=1 ; identity(3,3)=1 318 isym=indkk(ikpt2+(sppoldbl-1)*(isppol2-1)*nkpt2,2) 319 !DEBUG 320 !write(std_out,*)' wfconv : isym=',isym 321 !ENDDEBUG 322 itimrev=indkk(ikpt2+(sppoldbl-1)*(isppol2-1)*nkpt2,6) 323 if(isym/=0)then 324 symrel_conv(:,:)=symrel(:,:,isym) 325 call mati3inv(symrel_conv,symm) 326 shiftg(:)=indkk(ikpt2+(sppoldbl-1)*(isppol2-1)*nkpt2,3:5) 327 tnons_conv(:)=tnons(:,isym) 328 if(sum(tnons_conv(:)**2)>tol8)then 329 ! Need to compute phase factors associated with nonsymmorphic translations. 330 conv_tnons=1 331 ABI_ALLOCATE(phase3d,(2,npw1)) 332 ABI_ALLOCATE(phase1d,(2,(2*n1+1)+(2*n2+1)+(2*n3+1))) 333 ! Although the routine getph is originally written for 334 ! atomic phase factors, it does precisely what we want 335 atindx(1)=1 336 call getph(atindx,1,n1,n2,n3,phase1d,tnons_conv) 337 end if 338 if(nspinor1==2 .and. nspinor2==2)then 339 ! Compute rotation in spinor space 340 call getspinrot(rprimd2,spinrot,symrel_conv) 341 end if 342 else 343 shiftg(:)=0 344 symm(:,:)=identity(:,:) 345 spinrot(:)=zero 346 spinrot(1)=one 347 end if 348 if(itimrev/=0)then 349 symm(:,:)=-symm(:,:) 350 end if 351 352 !DEBUG 353 !write(std_out,'(a,i3,2x,3i3,2x,9i3)')' wfconv : isym,shiftg,symm=',isym,shiftg,symm 354 !write(std_out,*)' wfconv : ecut2_eff,ecut1=',ecut2_eff,ecut1 355 !write(std_out,*)' wfconv : istwf1_k,istwf2_k=',istwf1_k,istwf2_k 356 !write(std_out,*)' wfconv : kpoint1(:),kpoint2_sph(:)=',& 357 !& kpoint1(:),kpoint2_sph(:) 358 !write(std_out,*)' wfconv : nspinor1,nspinor2=',nspinor1,nspinor2 359 !ENDDEBUG 360 361 !if (mpi_enreg1%fft_option_lob==0) mpi_enreg1%fft_option_lob=1 362 !if (mpi_enreg2%fft_option_lob==0) mpi_enreg2%fft_option_lob=1 363 364 if (restart==2.and.(convert/=0.or.(nbd2/nspinor2>nbd1/nspinor1.and.formeig==0))) then 365 ! kg2 is needed both for FFT grid conversion and for envlop 366 ! Choose the center of the sphere : either gamma, or each k-point 367 kpoint2_sph(:)=0.0_dp 368 if(ceksp2==0)kpoint2_sph(:)=kptns2(:,ikpt2) 369 istwf2_k=istwfk2(ikpt2) 370 call kpgsph(ecut2_eff,exchn2n3d,gmet2,0,ikpt2,istwf2_k,kg2,kpoint2_sph,1,mpi_enreg2,mpw2,npw2) 371 end if 372 373 if(convert/=0)then 374 istwf10_k=0 375 if(ikpt10/=0)istwf10_k=istwfk1(ikpt10) 376 377 ! Only need G sphere if different from last time 378 if ( ikpt1/=ikpt10 .or. istwf1_k/=istwf10_k ) then 379 380 call kpgsph (ecut1,exchn2n3d,gmet1,0,ikpt1,istwf1_k,kg1,kpoint1,1,mpi_enreg1,mpw1,npw1) 381 if (debug>0) then 382 write(message, '(a,f8.3,a,a,3f8.5,a,a,i3,a,3(a,3es16.8,a),a,3i4,a,i5,a)' )& 383 & ' wfconv: called kpgsph with ecut1=',ecut1,ch10,& 384 & ' kpt1=',kptns1(1:3,ikpt1),ch10,& 385 & ' istwf1_k=',istwf1_k,ch10,& 386 & ' gmet1= ',gmet1(1:3,1),ch10,& 387 & ' ',gmet1(1:3,2),ch10,& 388 & ' ',gmet1(1:3,3),ch10,& 389 & ' ngfft=',ngfft_now(1:3),' giving npw1=',npw1,'.' 390 call wrtout(std_out,message,'PERS') 391 end if 392 ikpt10 = ikpt1 393 istwf10_k=istwf1_k 394 end if 395 396 if(conv_tnons==1)then 397 arg=two_pi*(kpoint1(1)*tnons_conv(1)+ kpoint1(2)*tnons_conv(2)+ kpoint1(3)*tnons_conv(3) ) 398 phktnons(1,1)=cos(arg) 399 phktnons(2,1)=sin(arg) 400 ! Convert 1D phase factors to 3D phase factors exp(i 2 pi (k+G).tnons ) 401 call ph1d3d(1,1,kg1,1,1,npw1,n1,n2,n3,phktnons,phase1d,phase3d) 402 end if 403 404 ABI_ALLOCATE(cfft,(2,n4,n5,n6)) 405 ABI_ALLOCATE(wavef1,(2,npw1)) 406 ABI_ALLOCATE(wavef2,(2,npw2)) 407 if(nspinor1==2 .and. nspinor2==2) then 408 ABI_ALLOCATE(wavefspinor,(2,2*npw2)) 409 end if 410 ABI_ALLOCATE(gbound1,(2*mgfft1+8,2)) 411 ABI_ALLOCATE(gbound2,(2*mgfft2+8,2)) 412 call sphereboundary(gbound1,istwf1_k,kg1,mgfft1,npw1) 413 call sphereboundary(gbound2,istwf2_k,kg2,mgfft2,npw2) 414 415 ! Take old wf from sphere->box, the new from box->sphere 416 ! One pays attention not to have a problem of erasing data when replacing 417 ! a small set of coefficient by a large set, or the reverse. 418 ! This is the reason of the use of order, _first and _last variables, 419 ! defined earlier. 420 nspinor_index=mpi_enreg1%me_spinor+1 421 do iband=iband_first,iband_last,order 422 do ispinor1=ispinor_first,ispinor_last,order 423 ispinor=ispinor1 424 if (mpi_enreg1%paral_spinor==1) then 425 if (ispinor1==nspinor_index) then 426 ispinor=1 427 else 428 if (nspinor1==2.and.nspinor2==2) wavefspinor(:,(ispinor1-1)*npw2+1:ispinor1*npw2)=zero 429 cycle 430 end if 431 end if 432 433 ! Copy input wf 434 i1=(ispinor-1)*npw1+(iband-1)*nspinor1_this_proc*npw1+icg1 435 wavef1(:,1:npw1)=cg1(:,i1+1:i1+npw1) 436 437 ! Make symmetry-induced conversion, if needed (translation part) 438 if(conv_tnons==1)then 439 !$OMP PARALLEL DO PRIVATE(ai,ar) 440 do ipw=1,npw1 441 ar=phase3d(1,ipw)*wavef1(1,ipw)-phase3d(2,ipw)*wavef1(2,ipw) 442 ai=phase3d(2,ipw)*wavef1(1,ipw)+phase3d(1,ipw)*wavef1(2,ipw) 443 wavef1(1,ipw)=ar 444 wavef1(2,ipw)=ai 445 end do 446 end if 447 448 ! Take into account time-reversal symmetry, if needed, in the scalar case 449 if(itimrev==1 .and. (nspinor1==1 .or. nspinor2==1))then 450 !$OMP PARALLEL DO 451 do ipw=1,npw1 452 wavef1(2,ipw)=-wavef1(2,ipw) 453 end do 454 end if 455 456 ! DEBUG 457 ! write(std_out,*)' wfconv : before sphere, isym,ispinor=',isym,ispinor 458 ! write(std_out,*)' no_shift,identity=',no_shift,identity 459 ! write(std_out,*)' shiftg,symm=',shiftg,symm 460 ! stop 461 ! This debugging sequence is an attempt to rotate spinors, 462 ! and works indeed for test13, when symmetry 9 is used ... 463 ! if(isym==9 .and. ispinor==1)then 464 ! write(std_out,*)' wfconv : gives a 120 degree rotation to first component' 465 ! do ipw=1,npw1 466 ! ar=- half*wavef1(1,ipw)-sqrt(three)*half*wavef1(2,ipw) 467 ! ai= sqrt(three)*half*wavef1(1,ipw)- half*wavef1(2,ipw) 468 ! wavef1(1,ipw)=ar 469 ! wavef1(2,ipw)=ai 470 ! end do 471 ! end if 472 ! ENDDEBUG 473 474 ! Convert wf, and also include the symmetry operation and shiftg. 475 call sphere(wavef1,1,npw1,cfft,n1,n2,n3,n4,n5,n6,kg1,istwf1_k,tobox,& 476 & mpi_enreg1%me_g0,no_shift,identity,one) 477 478 call sphere(wavef2,1,npw2,cfft,n1,n2,n3,n4,n5,n6,kg2,istwf2_k,tosph,& 479 & mpi_enreg2%me_g0,shiftg,symm,one) 480 481 if(nspinor2==1 )then 482 i2=(ispinor-1)*npw2+(iband-1)*nspinor2_this_proc*npw2+icg2 483 cg2(:,i2+1:i2+npw2)=wavef2(:,1:npw2) 484 else if(nspinor1==2.and.nspinor2==2)then 485 ! Will treat this case outside of the ispinor loop 486 i2=(ispinor1-1)*npw2 487 wavefspinor(:,i2+1:i2+npw2)=wavef2(:,1:npw2) 488 else if(nspinor1==1 .and. nspinor2==2)then 489 ! The number of bands is doubled, and the number of coefficients 490 ! is doubled also 491 if (mpi_enreg2%paral_spinor==0) then 492 i2=(iband-1)*nspinor2_this_proc*nspinor2_this_proc*npw2+icg2 493 cg2(:,i2+1:i2+npw2)=wavef2(:,1:npw2) 494 cg2(:,i2+npw2+1:i2+2*npw2)=zero 495 cg2(:,i2+2*npw2+1:i2+3*npw2)=zero 496 cg2(:,i2+3*npw2+1:i2+4*npw2)=wavef2(:,1:npw2) 497 else 498 i2=(iband-1)*nspinor2_this_proc*npw2+icg2 499 if (nspinor_index==1) then 500 cg2(:,i2+1:i2+npw2)=wavef2(:,1:npw2) 501 cg2(:,i2+npw2+1:i2+2*npw2)=zero 502 else 503 cg2(:,i2+1:i2+npw2)=zero 504 cg2(:,i2+npw2+1:i2+2*npw2)=wavef2(:,1:npw2) 505 end if 506 end if 507 end if 508 end do ! ispinor=ispinor_first,ispinor_last,order 509 510 if(nspinor1==2.and.nspinor2==2)then 511 ! Take care of possible parallelization over spinors 512 if (mpi_enreg2%paral_spinor==1) then 513 call xmpi_sum(wavefspinor,mpi_enreg2%comm_spinor,ierr) 514 end if 515 ! Take care of time-reversal symmetry, if needed 516 if(itimrev==1)then 517 ! Exchange spin-up and spin-down 518 ! Make complex conjugate of one component, 519 ! and change sign of other component 520 !$OMP PARALLEL DO PRIVATE(ipw,ar,ai) SHARED(wavefspinor,npw2) 521 do ipw=1,npw2 522 ! Here, change sign of real part 523 ar=-wavefspinor(1,ipw) 524 ai= wavefspinor(2,ipw) 525 wavefspinor(1,ipw)= wavefspinor(1,npw2+ipw) 526 ! Here, change sign of imaginary part 527 wavefspinor(2,ipw)=-wavefspinor(2,npw2+ipw) 528 wavefspinor(1,npw2+ipw)=ar 529 wavefspinor(2,npw2+ipw)=ai 530 end do 531 end if ! itimrev==1 532 533 ! Rotation in spinor space 534 !$OMP PARALLEL DEFAULT(PRIVATE) SHARED(npw2,spinrot,wavefspinor) 535 spinrots=spinrot(1) 536 spinrotx=spinrot(2) 537 spinroty=spinrot(3) 538 spinrotz=spinrot(4) 539 !$OMP DO 540 do ipw=1,npw2 541 ar=wavefspinor(1,ipw) 542 ai=wavefspinor(2,ipw) 543 br=wavefspinor(1,npw2+ipw) 544 bi=wavefspinor(2,npw2+ipw) 545 wavefspinor(1,ipw) = spinrots*ar-spinrotz*ai +spinroty*br-spinrotx*bi 546 wavefspinor(2,ipw) = spinrots*ai+spinrotz*ar +spinroty*bi+spinrotx*br 547 wavefspinor(1,npw2+ipw)=-spinroty*ar-spinrotx*ai +spinrots*br+spinrotz*bi 548 wavefspinor(2,npw2+ipw)=-spinroty*ai+spinrotx*ar +spinrots*bi-spinrotz*br 549 end do 550 !$OMP END DO 551 !$OMP END PARALLEL 552 553 ! Save wavefunction 554 i2=(iband-1)*nspinor2_this_proc*npw2+icg2 555 if (mpi_enreg2%paral_spinor==0) then 556 cg2(:,i2 +1:i2+ npw2)=wavefspinor(:,1:npw2) 557 cg2(:,i2+npw2+1:i2+2*npw2)=wavefspinor(:,npw2+1:2*npw2) 558 else 559 if (nspinor_index==1) then 560 cg2(:,i2+1:i2+npw2)=wavefspinor(:,1:npw2) 561 else 562 cg2(:,i2+1:i2+npw2)=wavefspinor(:,npw2+1:2*npw2) 563 end if 564 end if 565 end if ! nspinor1==2 .and. nspinor2==2 566 567 end do 568 569 ! Take care of copying eig and occ when nspinor increases or decreases 570 if(nspinor1==1.and.nspinor2==2)then 571 if(formeig==0)then 572 ! Note the reverse order, needed in case inplace=1 573 do iband=nbd1,1,-1 574 ! use eig_tmp to avoid bug on ifort10.1 x86_64 575 eig_tmp=eig_k1(iband) 576 eig_k2(2*iband-1:2*iband)=eig_tmp 577 ! occ_tmp=occ_k1(iband)*0.5_dp 578 ! occ_k2(2*iband-1:2*iband )=occ_tmp 579 end do 580 else 581 call wrtout(std_out,' wfconv: not yet coded, formeig=1!',"COLL") 582 end if 583 end if 584 if(nspinor1==2 .and. nspinor2==1)then 585 if(formeig==0)then 586 do iband=1,nbd1 587 ! use eig_tmp to avoid bug on ifort10.1 x86_64 588 eig_tmp=eig_k1(2*iband-1) 589 eig_k2(iband)=eig_tmp 590 ! occ_tmp=occ_k1(2*iband-1)*2.0_dp 591 ! occ_k2(iband)=occ_tmp 592 end do 593 else 594 call wrtout(std_out,' wfconv: not yet coded, formeig=1!',"COLL") 595 end if 596 end if 597 598 ABI_DEALLOCATE(cfft) 599 ABI_DEALLOCATE(gbound1) 600 ABI_DEALLOCATE(gbound2) 601 ABI_DEALLOCATE(wavef1) 602 ABI_DEALLOCATE(wavef2) 603 if(nspinor1==2 .and. nspinor2==2) then 604 ABI_DEALLOCATE(wavefspinor) 605 end if 606 607 else if(convert==0)then 608 609 if(inplace==0)then 610 ! Must copy cg, eig and occ if not in-place while convert==0 611 ! Note that npw1=npw2, nspinor1=nspinor2 612 cg2(:,1+icg2:npw1*nspinor1_this_proc*nbd1+icg2)=& 613 & cg1(:,1+icg1:npw1*nspinor1_this_proc*nbd1+icg1) 614 eig_k2(:)=eig_k1(:) 615 ! occ_k2(:)=occ_k1(:) 616 end if 617 618 end if ! End of if convert/=0 619 620 if(conv_tnons==1) then 621 ABI_DEALLOCATE(phase1d) 622 ABI_DEALLOCATE(phase3d) 623 end if 624 625 626 !If not enough bands, complete with random numbers or zeros 627 if(nbd2/nspinor2>nbd1/nspinor1)then 628 if(formeig==0)then 629 630 ! Ground state wf and eig case 631 eig_k2((nbd1/nspinor1)*nspinor2+1:nbd2)=huge(0.0_dp)/10.0_dp 632 occ_k2((nbd1/nspinor1)*nspinor2+1:nbd2)=0.0_dp 633 index=(nbd1/nspinor1)*nspinor2*npw2*nspinor2_this_proc 634 635 ! Initialisation of wavefunctions 636 ! One needs to initialize wfs in such a way to avoid symmetry traps, 637 ! and to avoid linear dependencies between wavefunctions 638 ! No need for a difference for different k points and/or spin-polarization 639 640 if (mpi_enreg1%paral_kgb == 1) then 641 npwtot=npw2 642 call timab(539,1,tsec) 643 call xmpi_sum(npwtot, mpi_enreg1%comm_bandfft, ierr) 644 call timab(539,2,tsec) 645 end if 646 647 do iband=(nbd1/nspinor1)*nspinor2+1,nbd2 648 do ispinor2=1,nspinor2_this_proc 649 ispinor=ispinor2;if (nspinor2_this_proc/=nspinor2) ispinor=mpi_enreg2%me_spinor+1 650 jsign=1;if (ispinor==2) jsign=-1 651 652 do ipw=1,npw2 653 index=index+1 654 ! Different seed for different planewave and band 655 ! DEBUG seq==par 656 ! if(.false.) then 657 ! ENDDEBUG seq==par 658 659 if ( mpi_enreg2%paral_kgb /= 1) then 660 seed=(iband-1)*npw2*nspinor2 + (ispinor-1)*npw2 + ipw 661 else 662 seed=jsign*(iband*(kg2(1,ipw)*npwtot*npwtot + kg2(2,ipw)*npwtot + kg2(3,ipw))) 663 end if 664 665 if(randalg == 0) then 666 ! For portability, use only integer numbers 667 ! The series of couples (fold1,fold2) is periodic with a period of 668 ! 3x5x7x11x13x17x19x23x29x31, that is, larger than 2**32, the largest integer*4 669 ! fold1 is between 0 and 34, fold2 is between 0 and 114. As sums of five 670 ! uniform random variables, their distribution is close to a gaussian 671 fold1=mod(seed,3)+mod(seed,5)+mod(seed,7)+mod(seed,11)+mod(seed,13) 672 fold2=mod(seed,17)+mod(seed,19)+mod(seed,23)+mod(seed,29)+mod(seed,31) 673 ! The gaussian distributions are folded, in order to be back to a uniform distribution 674 ! foldre is between 0 and 20, foldim is between 0 and 18 675 foldre=mod(fold1+fold2,21) 676 foldim=mod(3*fold1+2*fold2,19) 677 cg2(1,index+icg2)=dble(foldre) 678 cg2(2,index+icg2)=dble(foldim) 679 else 680 ! (AL) Simple linear congruential generator from 681 ! numerical recipes, modulo'ed and 64bit'ed to avoid 682 ! overflows (NAG doesn't like overflows, even though 683 ! they are perfectly legitimate here). Then, we get some 684 ! lowest order bits and sum them, as the previous 685 ! generator, to get quasi-normal numbers. 686 ! This is clearly suboptimal and might cause problems, 687 ! but at least it doesn't seem to create linear 688 ! dependencies and local minima like the previous one. 689 ! it's not trivial to generate good reproductible random 690 ! numbers in parallel. Patches welcome ! 691 ! Note a fun fortran fact : MOD simply ignores 64 bits integer 692 ! and casts them into 32bits, so we use MODULO. 693 fold1 = modulo(1664525_int64 * seed + 1013904223_int64, 2147483648_int64) 694 fold2 = modulo(1664525_int64 * fold1 + 1013904223_int64, 2147483648_int64) 695 fold1=modulo(fold1,3)+modulo(fold1,5)+modulo(fold1,7)+modulo(fold1,11)+modulo(fold1,13) 696 fold2=modulo(fold2,3)+modulo(fold2,5)+modulo(fold2,7)+modulo(fold2,11)+modulo(fold2,13) 697 cg2(1,index+icg2)=dble(fold1)/34-0.5 698 cg2(2,index+icg2)=dble(fold2)/34-0.5 699 end if 700 end do 701 end do 702 703 ! XG030513 : MPIWF need to impose cg to zero when at Gamma 704 ! Time-reversal symmetry for k=gamma impose zero imaginary part at G=0 705 ! XG : I do not know what happens for spin-orbit here : 706 if(istwf2_k==2 .and. mpi_enreg2%me_g0==1) then 707 cg2(2,1+(iband-1)*npw2*nspinor2_this_proc+icg2)=zero 708 end if 709 end do 710 711 ! Multiply with envelope function to reduce kinetic energy 712 icgmod=icg2+npw2*nspinor2_this_proc*(nbd1/nspinor1) 713 nbremn=nbd2-nbd1 714 call cg_envlop(cg2,ecut2,gmet2,icgmod,kg2,kpoint2_sph,mcg2,nbremn,npw2,nspinor2_this_proc) 715 716 if(ikpt2<=nkpt_max)then 717 write(message,'(3(a,i6))')' wfconv:',nbremn,' bands initialized randomly with npw=',npw2,', for ikpt=',ikpt2 718 call wrtout(std_out,message,'PERS') 719 end if 720 721 else if(formeig==1)then 722 723 ! For response function, put large numbers in the remaining of the 724 ! eigenvalue array (part of it was already filled in calling routine) 725 ! WARNING : Change of nspinor not yet coded 726 eig_k2(1+2*nbd1*nbd2 : 2*nbd2*nbd2)=huge(0.0_dp)/10.0_dp 727 ! Initialisation of wfs with 0 s 728 index=npw2*nbd1*nspinor2_this_proc 729 do iband=nbd1+1,nbd2 730 do ipw=1,npw2*nspinor2_this_proc 731 index=index+1 732 cg2(:,index+icg2)=zero 733 end do 734 end do 735 736 if(ikpt2<=nkpt_max)then 737 nbremn=nbd2-nbd1 738 write(message,'(a,i6,a,i7,a,i4)')' wfconv :',nbremn,' bands set=0 with npw=',npw2,', for ikpt=',ikpt2 739 call wrtout(std_out,message,'PERS') 740 end if 741 742 end if ! End of initialisation to 0 743 end if 744 745 !Orthogonalize GS wfs 746 !if (.False.) then 747 if (optorth==1.and.formeig==0.and.mpi_enreg2%paral_kgb/=1) then 748 ABI_ALLOCATE(dum,(2,0)) 749 ortalgo=0 !;ortalgo=3 750 call pw_orthon(icg2,0,istwf2_k,mcg2,0,npw2*nspinor2_this_proc,nbd2,ortalgo,dum,0,cg2,& 751 & mpi_enreg2%me_g0,mpi_enreg2%comm_bandspinorfft) 752 ABI_DEALLOCATE(dum) 753 end if 754 755 end subroutine wfconv