TABLE OF CONTENTS
ABINIT/pawmkaewf [ Functions ]
NAME
pawmkaewf
FUNCTION
Construct complete AE wave functions on the fine FFT grid adding onsite PAW corrections.
COPYRIGHT
Copyright (C) 2008-2018 ABINIT group (MG) This file is distributed under the terms of the GNU General Public License, see ~ABINIT/COPYING or http://www.gnu.org/copyleft/gpl.txt .
INPUTS
crystal<crystal_t>=Crystalline structure ebands<ebands_t>=Electronic energies dimcprj(natom)=array of dimensions of array cprj (not ordered) mband=maximum number of bands mcg=size of wave-functions array (cg) =mpw*nspinor*mband*mkmem*nsppol mcprj=size of projected wave-functions array (cprj) =nspinor*mband*mkmem*nsppol mkmem=number of k points treated by this node. mpi_atmtab(:)=--optional-- indexes of the atoms treated by current proc [comm_atom]= MPI communicator over atoms mpw=maximum dimensioned size of npw. my_natom=number of atoms treated by current processor natom=number of atoms in cell ntypat=number of types of atoms in the cell nkpt=Total number of k-points nsppol=1 for unpolarized, 2 for spin-polarized unks=unit number for G vectors. nband(nkpt*nsppol)=Number of bands for each k-point and spin. istwfk(nkpt)=Storage mode at each k-point. Pawfgrtab(natom) <type(pawfgrtab_type)> : data about the fine grid around each atom Pawrad(ntypat) <type(pawrad_type)> : radial mesh data for each type of atom Pawtab(ntypat) <type(pawtab_type)> : PAW functions around each type of atom Dtfil <type(datafiles_type)>=variables related to files cg(2,mcg)=planewave coefficients of wavefunctions. Cprj(natom,nspinor*mband*mkmem*nsppol)=<p_lmn|Cnk> coefficients for each WF |Cnk> and each |p_lmn> non-local projector npwarr(nkpt)=Number of plane waves at each k-point ngfftf(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft Note that ngfftf refers to the fine mesh. kg(3,mpw*mkmem)=reduced planewave coordinates Hdr<hdr_type>=the header of wf, den and pot files kpt(3,nkpt)=reduced coordinates of k points.
OUTPUT
ierr=Status error Main output is written on file (ETSF_IO file format).
NOTES
In PAW calculations, the pseudized wavefunction us represented on a relatively small plane wave basis set and is not normalized as it does not include the on-site PAW contributions which is described in terms of real spherical harmonics and radial functions. For post-processing and proper visualization, it is necessary to use the full electronic wave function, which is what this subroutine constructs. Specifically, it computes the pseudo part by doing an FFT from G- to r-space using the dense mesh defined by pawecutdg. The on-site PAW terms are also computed in real space inside each sphere and added to the pseudo part. Notice that this formula is expressed on the fine grid, and requires interpolating the PAW radial functions onto this grid, as well as calling initylmr in order to get the angular functions on the grid points.
PARENTS
outscfcv
CHILDREN
flush_unit,fourwf,free_my_atmtab,get_my_atmtab,nhatgrid paw_pwaves_lmn_free,paw_pwaves_lmn_init,pawcprj_alloc,pawcprj_free pawfgrtab_free,pawfgrtab_init,pawfgrtab_print,pawtab_get_lsize sphereboundary,wrap2_zero_one,wrtout,xcart2xred,xmpi_barrier,xmpi_max xmpi_sum
SOURCE
78 #if defined HAVE_CONFIG_H 79 #include "config.h" 80 #endif 81 82 #include "abi_common.h" 83 84 subroutine pawmkaewf(Dtset,crystal,ebands,my_natom,mpw,mband,mcg,mcprj,nkpt,mkmem,nsppol,nband,& 85 & istwfk,npwarr,kpt,ngfftf,kg,dimcprj,Pawfgrtab,Pawrad,Pawtab,& 86 & Hdr,Dtfil,cg,Cprj,MPI_enreg,ierr,pseudo_norms,set_k,set_band , & 87 & mpi_atmtab,comm_atom) ! Optional arguments 88 89 90 use defs_basis 91 use defs_datatypes 92 use defs_abitypes 93 use defs_wvltypes 94 use m_profiling_abi 95 use m_xmpi 96 use m_blas 97 use m_splines 98 use m_errors 99 use m_nctk 100 use m_hdr 101 #ifdef HAVE_NETCDF 102 use netcdf 103 #endif 104 105 use m_io_tools, only : flush_unit 106 use m_numeric_tools, only : wrap2_zero_one 107 use m_pptools, only : printxsf 108 use m_geometry, only : xcart2xred 109 use m_crystal, only : crystal_t 110 use m_crystal_io, only : crystal_ncwrite 111 use m_ebands, only : ebands_ncwrite 112 use m_pawrad, only : pawrad_type 113 use m_pawtab, only : pawtab_type, pawtab_get_lsize 114 use m_pawfgrtab, only : pawfgrtab_type, pawfgrtab_init, pawfgrtab_free, pawfgrtab_print 115 use m_pawcprj, only : pawcprj_type, pawcprj_alloc, pawcprj_get, pawcprj_free 116 use m_paw_pwaves_lmn, only : paw_pwaves_lmn_t, paw_pwaves_lmn_init, paw_pwaves_lmn_free 117 use m_paral_atom, only : get_my_atmtab, free_my_atmtab 118 119 !This section has been created automatically by the script Abilint (TD). 120 !Do not modify the following lines by hand. 121 #undef ABI_FUNC 122 #define ABI_FUNC 'pawmkaewf' 123 use interfaces_14_hidewrite 124 use interfaces_32_util 125 use interfaces_52_fft_mpi_noabirule 126 use interfaces_53_ffts 127 use interfaces_65_paw, except_this_one => pawmkaewf 128 !End of the abilint section 129 130 implicit none 131 132 !Arguments ------------------------------------ 133 !scalars 134 integer,intent(in) :: my_natom,mband,mcg,mcprj,mkmem,mpw,nsppol,nkpt 135 integer,intent(in),optional :: comm_atom,set_k,set_band 136 integer,intent(out) :: ierr 137 type(Datafiles_type),intent(in) :: Dtfil 138 type(MPI_type),intent(in) :: MPI_enreg 139 type(hdr_type),intent(inout) :: Hdr 140 type(dataset_type),intent(in) :: Dtset 141 type(crystal_t),intent(in) :: crystal 142 type(ebands_t),intent(in) :: ebands 143 !arrays 144 integer,intent(in) :: nband(nkpt*nsppol),istwfk(nkpt),npwarr(nkpt),dimcprj(crystal%natom) 145 integer,intent(in) :: ngfftf(18),kg(3,mpw*mkmem) 146 integer,optional,target,intent(in) :: mpi_atmtab(:) 147 real(dp),intent(in) :: cg(2,mcg) 148 real(dp),intent(in) :: kpt(3,nkpt) 149 real(dp),optional,intent(out) :: pseudo_norms(nsppol,nkpt,mband) 150 type(pawfgrtab_type),intent(in) :: Pawfgrtab(my_natom) 151 type(pawrad_type),intent(in) :: Pawrad(crystal%ntypat) 152 type(pawtab_type),intent(in) :: Pawtab(crystal%ntypat) 153 type(pawcprj_type),intent(in) :: Cprj(crystal%natom,mcprj) 154 155 !Local variables------------------------------- 156 !scalars 157 integer,parameter :: tim_fourwf0=0,tim_rwwf0=0,master=0 158 integer :: bdtot_index,iband,icg,mgfftf,paral_kgb 159 integer :: iatom,iatom_tot,ifgd,ifftsph,ifft,itypat,ispinor,ipw,ndat,ii,i1,i2,i3 160 integer :: jl,jm,jlmn,natom 161 integer :: max_nfgd,nfgd,ln_size,lmn_size,my_comm_atom,option 162 integer :: iorder_cprj,comm_cell,me_kpt,ibsp,ibg,isppol,ikpt,nband_k,cplex 163 integer :: n1,n2,n3,n4,n5,n6,ikg,npwout,istwf_k,npw_k 164 integer :: nfftot,nprocs 165 integer :: optcut,optgr0,optgr1,optgr2,optrad,start_band,start_kpt,stop_kpt,stop_band 166 logical :: my_atmtab_allocated,paral_atom 167 real(dp),parameter :: weight1=one 168 real(dp) :: phj,tphj,re_p,im_p,norm,norm_rerr,max_rerr,imur,reur,arg 169 character(len=500) :: msg 170 character(len=nctk_slen) :: shape_str 171 !arrays 172 integer,allocatable :: l_size_atm(:) 173 integer, pointer :: my_atmtab(:) 174 integer,allocatable :: gbound(:,:),kg_k(:,:) 175 real(dp) :: red(3),shift(3),rfft(3),kpoint(3),cp_fact(2) 176 real(dp),allocatable :: r0shift(:,:,:),phk_atm(:,:,:) 177 real(dp),allocatable :: buf_tmp(:,:,:),fofgin(:,:),fofgin_down(:,:),fofgout(:,:) 178 real(dp),allocatable :: denpot(:,:,:),fofr(:,:,:,:),fofr_down(:,:,:,:),phkr(:,:) 179 real(dp),allocatable :: ur_ae(:,:), ur_pw(:,:),ur_ae_onsite(:,:),ur_ps_onsite(:,:) 180 real(dp),allocatable :: ur_mask(:),dummy_1d(:),rsph_red(:,:),rsph_cart(:,:) 181 type(pawcprj_type),allocatable :: Cprj_k(:,:) 182 type(pawfgrtab_type) :: local_pawfgrtab(my_natom) 183 type(paw_pwaves_lmn_t),allocatable :: Paw_onsite(:) 184 #ifdef HAVE_NETCDF 185 integer :: fform,ncerr,ncid,ae_ncid,pw_ncid,aeons_ncid,psons_ncid 186 character(len=fnlen) :: fname 187 #endif 188 189 ! ************************************************************************ 190 191 DBG_ENTER("COLL") 192 193 !Init parallelism 194 comm_cell = MPI_enreg%comm_cell; nprocs = xmpi_comm_size(comm_cell) 195 me_kpt = MPI_enreg%me_kpt; paral_kgb=mpi_enreg%paral_kgb 196 197 !Compatibility tests 198 ABI_CHECK(mkmem/=0, "mkmem==0 not supported anymore!") 199 ABI_CHECK(MPI_enreg%paral_kgb == 0, "paral_kgb/=0 not coded") 200 ABI_CHECK(SIZE(dimcprj)>0, "dimcprj should be allocated") 201 ABI_CHECK(mpi_enreg%paral_spinor==0, "parallelisation over spinors not implemented") 202 ABI_CHECK(nprocs==1, "k spin parallelism not yet active") 203 ABI_CHECK(dtset%nspinor==1, "nspinor == 2 is buggy") 204 205 natom = crystal%natom 206 207 !Set up parallelism over atoms 208 paral_atom=(present(comm_atom).and.(my_natom/=natom)) 209 nullify(my_atmtab);if (present(mpi_atmtab)) my_atmtab => mpi_atmtab 210 my_comm_atom=xmpi_comm_self;if (present(comm_atom)) my_comm_atom=comm_atom 211 call get_my_atmtab(my_comm_atom,my_atmtab,my_atmtab_allocated,paral_atom,natom,my_natom_ref=my_natom) 212 213 !If collection of pseudo norms is enabled, make sure the array is initialised 214 if (present(pseudo_norms)) pseudo_norms = zero 215 216 !use a local copy of pawfgrtab to make sure we use the correction in the paw spheres 217 !the usual pawfgrtab uses r_shape which may not be the same as r_paw 218 if (my_natom>0) then 219 if (paral_atom) then 220 call pawtab_get_lsize(pawtab,l_size_atm,my_natom,Dtset%typat,mpi_atmtab=my_atmtab) 221 call pawfgrtab_init(local_pawfgrtab,Pawfgrtab(1)%cplex,l_size_atm,Dtset%nspden,Dtset%typat,& 222 & mpi_atmtab=my_atmtab,comm_atom=my_comm_atom) 223 else 224 call pawtab_get_lsize(pawtab,l_size_atm,my_natom,Dtset%typat) 225 call pawfgrtab_init(local_pawfgrtab,Pawfgrtab(1)%cplex,l_size_atm,Dtset%nspden,Dtset%typat) 226 end if 227 ABI_FREE(l_size_atm) 228 end if 229 optcut = 1 ! use rpaw to construct local_pawfgrtab 230 optgr0 = 0; optgr1 = 0; optgr2 = 0 ! dont need gY terms locally 231 optrad = 1 ! do store r-R 232 233 if (paral_atom) then 234 call nhatgrid(crystal%atindx1,crystal%gmet,my_natom,natom,crystal%nattyp,ngfftf,crystal%ntypat,& 235 & optcut,optgr0,optgr1,optgr2,optrad,local_pawfgrtab,pawtab,crystal%rprimd,Dtset%typat,crystal%ucvol,Hdr%xred,& 236 & comm_atom=my_comm_atom,mpi_atmtab=my_atmtab) 237 else 238 call nhatgrid(crystal%atindx1,crystal%gmet,my_natom,natom,crystal%nattyp,ngfftf,crystal%ntypat,& 239 & optcut,optgr0,optgr1,optgr2,optrad,local_pawfgrtab,pawtab,crystal%rprimd,Dtset%typat,crystal%ucvol,Hdr%xred) 240 end if 241 !now local_pawfgrtab is ready to use 242 243 max_nfgd=MAXVAL(local_pawfgrtab(:)%nfgd) ! MAX no. of points in the fine grid for this PAW sphere 244 ABI_MALLOC(r0shift,(3,max_nfgd,my_natom)) 245 ABI_MALLOC(phk_atm,(2,max_nfgd,my_natom)) 246 247 do iatom=1,my_natom 248 iatom_tot=iatom;if (paral_atom) iatom_tot=my_atmtab(iatom) 249 250 nfgd=local_pawfgrtab(iatom)%nfgd ! no. of points in the fine grid for this PAW sphere 251 ABI_MALLOC(rsph_red,(3,nfgd)) 252 ABI_MALLOC(rsph_cart,(3,nfgd)) 253 do ifgd=1,nfgd 254 rsph_cart(:,ifgd) = local_pawfgrtab(iatom)%rfgd(:,ifgd) + crystal%xcart(:,iatom_tot) 255 end do 256 call xcart2xred(nfgd,crystal%rprimd,rsph_cart,rsph_red) ! we work in reduced coordinates. 257 do ifgd=1,nfgd 258 call wrap2_zero_one(rsph_red(1,ifgd),red(1),shift(1)) ! num = red + shift 259 call wrap2_zero_one(rsph_red(2,ifgd),red(2),shift(2)) 260 call wrap2_zero_one(rsph_red(3,ifgd),red(3),shift(3)) 261 r0shift(:,ifgd,iatom) = shift 262 !if (ANY( ABS(shift) > tol12)) then 263 ! MSG_WARNING("rmR_red is outside the first unit cell.") 264 ! write(std_out,*)rsph_red(:,ifgd),shift 265 !end if 266 end do 267 ABI_FREE(rsph_red) 268 ABI_FREE(rsph_cart) 269 end do 270 271 if (.not.paral_atom .and. my_natom>0) then 272 call pawfgrtab_print(local_pawfgrtab,natom=natom,unit=std_out,& 273 & prtvol=Dtset%prtvol,mode_paral="COLL") 274 end if 275 276 ierr=0 277 #ifndef HAVE_NETCDF 278 ierr = -1 279 write(msg,'(3a)')& 280 & "netcdf support must be enabled in order to output AE PAW wavefunction. ",ch10,& 281 & "No output will be produced, use --enable-netcdf at configure-time. " 282 MSG_WARNING(msg) 283 return 284 !These statements are necessary to avoid the compiler complain about unused variables: 285 ii=Dtset%usepaw;ii=Dtfil%unpaw;ii=Hdr%usepaw 286 #endif 287 288 !FIXME check ordering in cprj and Eventually in external file 289 !why is iorder_cprj not stored in the file for crosschecking purpose? 290 !Here Im assuming cprj are not ordered! 291 iorder_cprj=0 292 293 !n4,n5,n6 are FFT dimensions, modified to avoid cache trashing 294 n1=ngfftf(1); n2=ngfftf(2); n3=ngfftf(3) 295 n4=ngfftf(4); n5=ngfftf(5); n6=ngfftf(6) 296 nfftot=PRODUCT(ngfftf(1:3)) 297 mgfftf=MAXVAL(ngfftf(1:3)) 298 299 ABI_MALLOC(phkr,(2,nfftot)) 300 ABI_MALLOC(gbound,(2*mgfftf+8,2)) 301 302 #ifdef HAVE_NETCDF 303 !=== Initialize ETSF_IO files === 304 ! FIXME: nspinor == 2 is buggy 305 306 fname = trim(dtfil%filnam_ds(4))//'_PAWAVES.nc' 307 write(msg,'(2a)')' Opening file for AE PAW wave functions: ',trim(fname) 308 call wrtout(std_out,msg,'PERS') 309 call wrtout(ab_out,msg,'PERS') 310 311 if (xmpi_comm_rank(comm_cell) == master) then 312 NCF_CHECK(nctk_open_create(ncid, fname, xmpi_comm_self)) 313 314 fform=602 315 NCF_CHECK(hdr_ncwrite(hdr, ncid, fform, nc_define=.True.)) 316 317 ! Define wavefunctions in real space on the dense FFT mesh 318 ! Fortran layout: 319 !real_space_wavefunctions: double 8d array with shape: 320 ! [real_or_complex_wavefunctions] 321 ! [number_of_grid_points_vector1][number_of_grid_points_vector2][number_of_grid_points_vector3] 322 ! [number_of_spinor_components] 323 ! [max_number_of_states][number_of_kpoints][number_of_spins] 324 325 ncerr = nctk_def_dims(ncid, [ & 326 nctkdim_t("real_or_complex_wavefunctions", 2), & 327 nctkdim_t("number_of_grid_points_vector1", n1), & 328 nctkdim_t("number_of_grid_points_vector2", n2), & 329 nctkdim_t("number_of_grid_points_vector3", n3) & 330 ], defmode=.True.) 331 NCF_CHECK(ncerr) 332 333 shape_str = "real_or_complex_wavefunctions, & 334 & number_of_grid_points_vector1, number_of_grid_points_vector2, number_of_grid_points_vector3, & 335 & number_of_spinor_components, & 336 & max_number_of_states, number_of_kpoints, number_of_spins" 337 338 ! Define wavefunctions in real space. 339 ncerr = nctk_def_arrays(ncid, [& 340 nctkarr_t('ur_ae', "dp", shape_str),& 341 nctkarr_t('ur_pw', "dp", shape_str),& 342 nctkarr_t('ur_ae_onsite', "dp", shape_str),& 343 nctkarr_t('ur_ps_onsite', "dp", shape_str) & 344 ], defmode=.True.) 345 NCF_CHECK(ncerr) 346 347 ! Complete the geometry information. 348 NCF_CHECK(crystal_ncwrite(crystal, ncid)) 349 NCF_CHECK(ebands_ncwrite(ebands, ncid)) 350 351 ! Close the file. 352 NCF_CHECK(nf90_close(ncid)) 353 end if 354 355 call xmpi_barrier(comm_cell) 356 357 ! Reopen the file in parallel inside comm_cell 358 ! Note that we use individual IO thus there's no need to handle idle processes 359 ! if paral_kgb == 0 and nprocs > nkpt * nsppol 360 NCF_CHECK(nctk_open_modify(ncid, fname, comm_cell)) 361 ae_ncid = nctk_idname(ncid, "ur_ae") 362 pw_ncid = nctk_idname(ncid, "ur_pw") 363 aeons_ncid = nctk_idname(ncid, "ur_ae_onsite") 364 psons_ncid = nctk_idname(ncid, "ur_ps_onsite") 365 366 NCF_CHECK(nctk_set_datamode(ncid)) 367 #endif 368 369 !Init structure storing phi_{nlm} and tphi_(nlm} on the dense FFT points located in the PAW spheres. 370 ABI_DT_MALLOC(Paw_onsite,(natom)) 371 if (paral_atom) then 372 call paw_pwaves_lmn_init(Paw_onsite,my_natom,natom,crystal%ntypat,crystal%rprimd,crystal%xcart,& 373 Pawtab,Pawrad,local_pawfgrtab, comm_atom=my_comm_atom,mpi_atmtab=my_atmtab) 374 else 375 call paw_pwaves_lmn_init(Paw_onsite,my_natom,natom,crystal%ntypat,crystal%rprimd,crystal%xcart,& 376 Pawtab,Pawrad,local_pawfgrtab) 377 end if 378 379 bdtot_index=0; icg=0; ibg=0; norm_rerr=smallest_real 380 381 ! === Loop over spin === 382 do isppol=1,nsppol 383 ikg=0; start_kpt=1; stop_kpt=nkpt 384 385 ! Check if k-point was specified (only serial) 386 if (present(set_k) .and. nprocs==1) then 387 if (set_k/=0) then 388 start_kpt = set_k 389 stop_kpt = set_k 390 !MSG_ERROR("set_k") 391 end if 392 end if 393 394 ! === Loop over k points === 395 do ikpt=start_kpt,stop_kpt 396 kpoint = kpt(:,ikpt) 397 nband_k = nband(ikpt+(isppol-1)*nkpt) 398 npw_k = npwarr(ikpt) 399 istwf_k = istwfk(ikpt) 400 401 if (proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,isppol,me_kpt)) then 402 bdtot_index=bdtot_index+nband_k 403 !MSG_ERROR("cycle in seq!") 404 cycle 405 end if 406 407 do i3=0,n3-1 408 rfft(3)=DBLE(i3)/n3 409 do i2=0,n2-1 410 rfft(2)=DBLE(i2)/n2 411 do i1=0,n1-1 412 rfft(1)=DBLE(i1)/n1 413 ifft = 1 +i1 +i2*n1 +i3*n1*n2 414 phkr(1,ifft) = COS(two_pi*dot_product(kpoint,rfft)) 415 phkr(2,ifft) = SIN(two_pi*dot_product(kpoint,rfft)) 416 end do 417 end do 418 end do 419 ! phkr(1,:)=one; phkr(2,:)=zero 420 421 ! Calculate the phase for the onsite PAW contributions. 422 do iatom=1,my_natom 423 nfgd=local_pawfgrtab(iatom)%nfgd ! no. of points in the fine grid for this PAW sphere 424 do ifgd=1,nfgd 425 arg = -two_pi* dot_product(r0shift(:,ifgd,iatom),kpoint) 426 phk_atm(1,ifgd,iatom) = COS(arg) 427 phk_atm(2,ifgd,iatom) = SIN(arg) 428 end do 429 end do 430 431 ABI_DT_MALLOC(Cprj_k,(natom,dtset%nspinor*nband_k)) 432 call pawcprj_alloc(Cprj_k,0,dimcprj) 433 434 ! Extract cprj for this k-point. 435 ibsp=0 436 do iband=1,nband_k 437 do ispinor=1,dtset%nspinor 438 ibsp=ibsp+1 439 do iatom=1,natom 440 Cprj_k(iatom,ibsp)%cp(:,:)=Cprj(iatom,ibsp+ibg)%cp(:,:) 441 end do 442 end do 443 end do 444 445 ABI_MALLOC(kg_k,(3,npw_k)) 446 447 ! Extract G-vectors. 448 kg_k(:,1:npw_k)=kg(:,1+ikg:npw_k+ikg) 449 call sphereboundary(gbound,istwf_k,kg_k,mgfftf,npw_k) 450 451 ! If a single band is requested, neuter the loop (only serial) 452 start_band = 1; stop_band = nband_k 453 if (present(set_band).AND.nprocs==1) then 454 if (set_band/=0) then 455 start_band = set_band 456 stop_band = set_band 457 !MSG_ERROR("set_band") 458 end if 459 end if 460 461 ! Loop over bands. 462 do iband=start_band,stop_band 463 464 ! Fourier transform on the real fft box of the smooth part. 465 ndat=Dtset%nspinor 466 ABI_MALLOC(fofgin,(2,npw_k*ndat)) 467 ABI_MALLOC(fofr,(2,n4,n5,n6*ndat)) 468 469 do ipw=1,npw_k*dtset%nspinor 470 fofgin(:,ipw)=cg(:,ipw+(iband-1)*npw_k*dtset%nspinor+icg) 471 end do 472 473 ! Complex can be set to 0 with this option(0) of fourwf 474 option=0; cplex=0; npwout=1 475 ABI_MALLOC(denpot,(cplex*n4,n5,n6)) 476 ABI_MALLOC(fofgout,(2,npwout*ndat)) 477 478 call fourwf(cplex,denpot,fofgin(:,1:npw_k),fofgout,fofr(:,:,:,1:n6),gbound,gbound,istwf_k,kg_k,kg_k,& 479 & mgfftf,MPI_enreg,1,ngfftf,npw_k,npwout,n4,n5,n6,option,paral_kgb,tim_fourwf0,weight1,weight1,& 480 & use_gpu_cuda=Dtset%use_gpu_cuda) 481 482 ! Here I do not know if fourwf works in the case of spinors, 483 ! It seems that not all fftalg option support ndata! should check! 484 ! Do not forget to declare real(dp)::fofgin_down(:,:) to use the following statements 485 if (Dtset%nspinor==2) then 486 ABI_MALLOC(fofgin_down,(2,npw_k)) 487 ABI_MALLOC(fofr_down,(2,n4,n5,n6)) 488 fofgin_down(:,:)=fofgin(:,1+npw_k:2*npw_k) 489 ! Complex can be set to 0 with this option(0) of fourwf 490 ! cplex=1; option=1; npwout=1; ndat=1 491 ! NOTE: fofr_down can NOT be replaced by fofr(:,:,:,n6+1:2*n6), or else 492 ! the data in fofr(:,:,:,1:n6) will be the same with fofr(:,:,:,n6+1:2*n6) 493 call fourwf(cplex,denpot,fofgin_down,fofgout,fofr_down,gbound,gbound,istwf_k,kg_k,kg_k,& 494 & mgfftf,MPI_enreg,1,ngfftf,npw_k,npwout,n4,n5,n6,option,paral_kgb,tim_fourwf0,weight1,weight1) 495 ABI_FREE(fofgin_down) 496 end if 497 498 ABI_MALLOC(ur_ae,(2,n1*n2*n3*ndat)) 499 ABI_MALLOC(ur_ae_onsite,(2,n1*n2*n3)) 500 ABI_MALLOC(ur_ps_onsite,(2,n1*n2*n3)) 501 ABI_MALLOC(ur_pw,(2,n1*n2*n3*ndat)) 502 ABI_MALLOC(ur_mask,(n1*n2*n3)) 503 504 ur_ae=zero;ur_ae_onsite=zero;ur_ps_onsite=zero;ur_pw=zero;ur_mask=zero 505 506 ! * Add phase e^{ikr} since it is contained in cprj. 507 do i3=1,n3 508 do i2=1,n2 509 do i1=1,n1 510 ii = i1 + n1*(i2-1)+ n1*n2*(i3-1) 511 ur_pw(:,ii)=fofr(:,i1,i2,i3) ! Save pw part separately without the phase. 512 ur_ae(1,ii)= fofr(1,i1,i2,i3) * phkr(1,ii) - fofr(2,i1,i2,i3) * phkr(2,ii) 513 ur_ae(2,ii)= fofr(1,i1,i2,i3) * phkr(2,ii) + fofr(2,i1,i2,i3) * phkr(1,ii) 514 if(Dtset%nspinor==2) then 515 ur_pw(:,ii+n1*n2*n3)=fofr_down(:,i1,i2,i3) ! Save pw part separately without the phase. 516 ur_ae(1,ii+n1*n2*n3)= fofr_down(1,i1,i2,i3) * phkr(1,ii) - fofr_down(2,i1,i2,i3) * phkr(2,ii) 517 ur_ae(2,ii+n1*n2*n3)= fofr_down(1,i1,i2,i3) * phkr(2,ii) + fofr_down(2,i1,i2,i3) * phkr(1,ii) 518 end if 519 end do 520 end do 521 end do 522 ABI_FREE(fofr) 523 524 if(Dtset%nspinor==2) then 525 ABI_FREE(fofr_down) 526 end if 527 528 ! === Add onsite term on the augmented FFT mesh === 529 do iatom=1,my_natom 530 itypat =local_pawfgrtab(iatom)%itypat 531 lmn_size=Pawtab(itypat)%lmn_size 532 ln_size =Pawtab(itypat)%basis_size ! no. of nl elements in PAW basis 533 nfgd =local_pawfgrtab(iatom)%nfgd ! no. of points in the fine grid for this PAW sphere 534 535 ibsp=(iband-1)*dtset%nspinor 536 do ispinor=1,dtset%nspinor 537 ibsp=ibsp+1 538 do jlmn=1,lmn_size 539 jl=Pawtab(itypat)%indlmn(1,jlmn) 540 jm=Pawtab(itypat)%indlmn(2,jlmn) 541 cp_fact(1) = Cprj_k(iatom,ibsp)%cp(1,jlmn) *sqrt(crystal%ucvol) ! Magic factor 542 cp_fact(2) = Cprj_k(iatom,ibsp)%cp(2,jlmn) *sqrt(crystal%ucvol) 543 544 do ifgd=1,nfgd ! loop over fine grid points in current PAW sphere. 545 ifftsph = local_pawfgrtab(iatom)%ifftsph(ifgd) ! index of the point on the grid 546 phj = Paw_onsite(iatom)% phi(ifgd,jlmn) 547 tphj = Paw_onsite(iatom)%tphi(ifgd,jlmn) 548 ! old code 549 !re_p = cp_fact(1); im_p = cp_fact(2) 550 ! apply the phase 551 re_p = cp_fact(1) * phk_atm(1,ifgd,iatom) - cp_fact(2) * phk_atm(2,ifgd,iatom) 552 im_p = cp_fact(1) * phk_atm(2,ifgd,iatom) + cp_fact(2) * phk_atm(1,ifgd,iatom) 553 554 ur_ae(1,ifftsph+(ispinor-1)*nfftot) = ur_ae(1,ifftsph+(ispinor-1)*nfftot) + re_p * (phj-tphj) 555 ur_ae(2,ifftsph+(ispinor-1)*nfftot) = ur_ae(2,ifftsph+(ispinor-1)*nfftot) + im_p * (phj-tphj) 556 ur_ae_onsite(1,ifftsph) = ur_ae_onsite(1,ifftsph) + re_p * phj 557 ur_ae_onsite(2,ifftsph) = ur_ae_onsite(2,ifftsph) + im_p * phj 558 ur_ps_onsite(1,ifftsph) = ur_ps_onsite(1,ifftsph) + re_p * tphj 559 ur_ps_onsite(2,ifftsph) = ur_ps_onsite(2,ifftsph) + im_p * tphj 560 ur_mask(ifftsph) = one 561 end do 562 563 end do !jlmn 564 end do !ispinor 565 end do !iatom 566 567 if (paral_atom) then 568 ABI_MALLOC(buf_tmp,(2,n1*n2*n3,3)) 569 buf_tmp(:,:,1) = ur_ae 570 buf_tmp(:,:,2) = ur_ae_onsite 571 buf_tmp(:,:,3) = ur_ps_onsite 572 call xmpi_sum(buf_tmp,my_comm_atom,ierr) 573 ur_ae = buf_tmp(:,:,1) 574 ur_ae_onsite= buf_tmp(:,:,2) 575 ur_ps_onsite= buf_tmp(:,:,3) 576 ABI_FREE(buf_tmp) 577 end if 578 579 ! * Remove the phase e^{ikr}, we store u(r). 580 do i3=1,n3 581 do i2=1,n2 582 do i1=1,n1 583 ii = i1 + n1*(i2-1)+ n1*n2*(i3-1) 584 reur=ur_ae(1,ii) 585 imur=ur_ae(2,ii) 586 ur_ae(1,ii)= reur * phkr(1,ii) + imur * phkr(2,ii) 587 ur_ae(2,ii)= -reur * phkr(2,ii) + imur * phkr(1,ii) 588 if(Dtset%nspinor==2) then 589 reur=ur_ae(1,ii+nfftot) ! Important! 590 imur=ur_ae(2,ii+nfftot) 591 ur_ae(1,ii+nfftot)= reur * phkr(1,ii) + imur * phkr(2,ii) 592 ur_ae(2,ii+nfftot)= -reur * phkr(2,ii) + imur * phkr(1,ii) 593 end if 594 reur=ur_ae_onsite(1,ii) 595 imur=ur_ae_onsite(2,ii) 596 ur_ae_onsite(1,ii)= reur * phkr(1,ii) + imur * phkr(2,ii) 597 ur_ae_onsite(2,ii)= -reur * phkr(2,ii) + imur * phkr(1,ii) 598 reur=ur_ps_onsite(1,ii) 599 imur=ur_ps_onsite(2,ii) 600 ur_ps_onsite(1,ii)= reur * phkr(1,ii) + imur * phkr(2,ii) 601 ur_ps_onsite(2,ii)= -reur * phkr(2,ii) + imur * phkr(1,ii) 602 end do 603 end do 604 end do 605 606 norm=zero 607 do ii=1,npw_k*Dtset%nspinor 608 norm=norm+fofgin(1,ii)**2+fofgin(2,ii)**2 609 end do 610 write(std_out,'(a,2i5,f22.16)',advance='no') 'ikpt,iband, norm (G,PSWF)=',ikpt,iband,norm 611 norm=zero 612 do ifft=1,nfftot*Dtset%nspinor 613 norm = norm + ur_ae(1,ifft)**2+ur_ae(2,ifft)**2 614 end do 615 norm=norm/nfftot 616 norm_rerr = MAX((ABS(norm-one))*100,norm_rerr) 617 write(std_out,*)"norm (R,AEWF)= ",norm 618 call flush_unit(std_out) 619 620 ! MS: Various testing and debugging options 621 if (.TRUE..and.nprocs==1) then 622 if (present(pseudo_norms)) then 623 ! Check the supposedly zero overlap |\tilde{Psi_n}-\tilde{Psi_n^1}|^2 624 ABI_MALLOC(dummy_1d,(n1*n2*n3)) 625 dummy_1d = zero 626 norm = zero 627 do ifft = 1, nfftot 628 dummy_1d(ifft) = ((ur_pw(1,ifft)-ur_ps_onsite(1,ifft))**2 & 629 & + (ur_pw(2,ifft)-ur_ps_onsite(2,ifft))**2) * ur_mask(ifft) 630 norm = norm + dummy_1d(ifft) 631 end do 632 norm = norm / nfftot 633 pseudo_norms(isppol,ikpt,iband) = norm 634 ABI_FREE(dummy_1d) 635 end if 636 637 else 638 write(msg,'(5a)')& 639 & "The option to print PAW all-electron wavefunctions is on, but execution ",ch10,& 640 & "is in parallel on two or more processors. XcrysDen files with individual con-",ch10,& 641 & "tributions will not be written. In order to enable this you must run in serial." 642 MSG_WARNING(msg) 643 end if ! Check if serial run 644 645 #ifdef HAVE_NETCDF 646 ncerr = nf90_put_var(ncid, ae_ncid, ur_ae, & 647 & start=[1,1,1,1,1,iband,ikpt,isppol], count=[2,n1,n2,n3,1,1,1,1]) 648 NCF_CHECK(ncerr) 649 650 ncerr = nf90_put_var(ncid, pw_ncid, ur_pw, & 651 & start=[1,1,1,1,1,iband,ikpt,isppol], count=[2,n1,n2,n3,1,1,1,1]) 652 NCF_CHECK(ncerr) 653 654 ncerr = nf90_put_var(ncid, aeons_ncid, ur_ae_onsite, & 655 & start=[1,1,1,1,1,iband,ikpt,isppol], count=[2,n1,n2,n3,1,1,1,1]) 656 NCF_CHECK(ncerr) 657 658 ncerr = nf90_put_var(ncid, psons_ncid, ur_ps_onsite, & 659 & start=[1,1,1,1,1,iband,ikpt,isppol], count=[2,n1,n2,n3,1,1,1,1]) 660 NCF_CHECK(ncerr) 661 #endif 662 663 ABI_FREE(ur_ae) 664 ABI_FREE(ur_ae_onsite) 665 ABI_FREE(ur_ps_onsite) 666 ABI_FREE(ur_pw) 667 ABI_FREE(ur_mask) 668 669 ABI_FREE(fofgin) 670 ABI_FREE(fofgout) 671 ABI_FREE(denpot) 672 end do !nband_k 673 674 bdtot_index=bdtot_index+nband_k 675 676 if (mkmem/=0) then 677 ibg=ibg+dtset%nspinor*nband_k 678 icg=icg+npw_k*dtset%nspinor*nband_k 679 ikg=ikg+npw_k 680 end if 681 682 ABI_FREE(kg_k) 683 684 call pawcprj_free(Cprj_k) 685 ABI_DT_FREE(Cprj_k) 686 687 end do !ikpt 688 end do !nsppol 689 690 ABI_FREE(phkr) 691 ABI_FREE(gbound) 692 693 ! Free augmentation waves. 694 call paw_pwaves_lmn_free(Paw_onsite) 695 ABI_DT_FREE(Paw_onsite) 696 697 ! Maximum relative error over CPUs. 698 call xmpi_max(norm_rerr,max_rerr,comm_cell,ierr) 699 write(std_out,*)"max_rerr=",max_rerr 700 701 if (max_rerr > ten) then 702 write(msg,'(7a)')& 703 & "Inaccuracy on the normalization of the wave funtions exceeds 10%. ",ch10,& 704 & "Likely due to the use of a too coarse FFT mesh or unconverged wavefunctions. ",ch10,& 705 & "Numerical values inside the augmentation regions might be inaccurate. ",ch10,& 706 & "Action: increase pawecutdg in your input file. " 707 MSG_COMMENT(msg) 708 end if 709 710 ABI_FREE(r0shift) 711 ABI_FREE(phk_atm) 712 call pawfgrtab_free(local_pawfgrtab) 713 714 !Destroy atom table used for parallelism 715 call free_my_atmtab(my_atmtab,my_atmtab_allocated) 716 717 DBG_EXIT("COLL") 718 719 end subroutine pawmkaewf