TABLE OF CONTENTS
ABINIT/denfgr [ Functions ]
NAME
denfgr
FUNCTION
Construct complete electron density on fine grid, by removing nhat and adding PAW corrections
COPYRIGHT
Copyright (C) 2005-2018 ABINIT group (JWZ) This file is distributed under the terms of the GNU General Public License, see ~ABINIT/COPYING or http://www.gnu.org/copyleft/gpl.txt .
INPUTS
atindx1(natom)=index table for atoms, inverse of atindx (see gstate.f) gmet(3,3)=reciprocal space metric tensor in bohr**-2. mpi_atmtab(:)=--optional-- indexes of the atoms treated by current proc comm_atom=--optional-- MPI communicator over atoms my_natom=number of atoms treated by current processor natom= number of atoms in cell nattyp(ntypat)= # atoms of each type. ngfft(18)=contain all needed information about 3D FFT (see NOTES at beginning of scfcv) nhat(pawfgr%nfft,nspden)= compensation charge density used in PAW nspinor=Number of spinor components nsppol=Number of independent spin components. nspden= number of spin densities ntypat= number of types of atoms in the cell pawfgr <type(pawfgr_type)>= data about the fine grid pawrad(ntypat) <type(pawrad_type)>= radial mesh data for each type of atom pawrhoij(natom) <type(pawrhoij_type)>= rho_ij data for each atom pawtab(ntypat) <type(pawtab_type)>= PAW functions around each type of atom rhor(pawfgr%nfft,nspden)= input density ($\tilde{n}+\hat{n}$ in PAW case) rprimd(3,3)=dimensional primitive translations for real space (bohr) typat(natom)= list of atom types ucvol=unit cell volume (bohr**3) xred(3,natom)=reduced dimensionless atomic coordinates
OUTPUT
rhor_paw(pawfgr%nfft,nspden)= full electron density on the fine grid
NOTES
In PAW calculations, the valence density present in rhor includes the
compensation charge density $\hat{n}$, and also doesn't include the on-site
PAW contributions. For post-processing and proper visualization it is necessary
to use the full electronic density, which is what this subroutine constructs.
Specifically, it removes $\hat{n}$ from rhor, and also computes the on-site PAW
terms. This is nothing other than the proper PAW treatment of the density
operator $|\mathbf{r}\rangle\langle\mathbf{r}|$, and yields the formula
$$\tilde{n}+\sum_{ij}\rho_ij\left[\varphi_i(\mathbf{r})\varphi_j(\mathbf{r})-
\tilde{\varphi}_i(\mathbf{r})\tilde{\varphi}_j(\mathbf{r})\right]$$
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
bethe_salpeter,outscfcv,sigma
CHILDREN
free_my_atmtab,get_my_atmtab,initylmr,nhatgrid,pawfgrtab_free pawfgrtab_init,pawrad_deducer0,pawtab_get_lsize,printxsf,sort_dp,spline splint,wrtout,xmpi_barrier,xmpi_sum,xred2xcart
SOURCE
67 #if defined HAVE_CONFIG_H 68 #include "config.h" 69 #endif 70 71 #include "abi_common.h" 72 73 subroutine denfgr(atindx1,gmet,spaceComm_in,my_natom,natom,nattyp,ngfft,nhat,nspinor,nsppol,nspden,ntypat, & 74 & pawfgr,pawrad,pawrhoij,pawtab,prtvol,rhor,rhor_paw,rhor_n_one,rhor_nt_one,rprimd,typat,ucvol,xred,& 75 & abs_n_tilde_nt_diff,znucl,mpi_atmtab,comm_atom) ! Optional arguments 76 77 use defs_basis 78 use m_profiling_abi 79 use m_errors 80 use m_xmpi 81 use m_splines 82 use m_sort 83 84 use m_io_tools, only : open_file 85 use m_geometry, only : xred2xcart 86 use m_pptools, only : printxsf 87 use m_pawrad, only : pawrad_type, pawrad_deducer0 88 use m_pawtab, only : pawtab_type,pawtab_get_lsize 89 use m_pawfgrtab, only : pawfgrtab_type, pawfgrtab_init, pawfgrtab_free 90 use m_pawrhoij, only : pawrhoij_type 91 use m_pawfgr, only : pawfgr_type 92 use m_paw_sphharm, only : initylmr 93 use m_paral_atom, only : get_my_atmtab, free_my_atmtab 94 95 !This section has been created automatically by the script Abilint (TD). 96 !Do not modify the following lines by hand. 97 #undef ABI_FUNC 98 #define ABI_FUNC 'denfgr' 99 use interfaces_14_hidewrite 100 use interfaces_65_paw, except_this_one => denfgr 101 !End of the abilint section 102 103 implicit none 104 105 !Arguments ------------------------------------ 106 !scalars 107 integer,intent(in) :: my_natom,natom,nspden,ntypat,prtvol,nsppol,nspinor 108 integer,optional,intent(in) :: comm_atom 109 real(dp),intent(in) :: ucvol 110 type(pawfgr_type),intent(in) :: pawfgr 111 !arrays 112 integer,intent(in) :: spaceComm_in 113 integer,intent(in) :: atindx1(natom),nattyp(ntypat),ngfft(18),typat(natom) 114 integer,optional,target,intent(in) :: mpi_atmtab(:) 115 real(dp),intent(in) :: gmet(3,3),nhat(pawfgr%nfft,nspden) 116 real(dp),intent(in) :: rhor(pawfgr%nfft,nspden),rprimd(3,3) 117 real(dp),intent(inout) :: xred(3,natom) 118 real(dp),intent(out) :: rhor_paw(pawfgr%nfft,nspden) 119 real(dp),intent(out) :: rhor_n_one(pawfgr%nfft,nspden) 120 real(dp),intent(out) :: rhor_nt_one(pawfgr%nfft,nspden) 121 real(dp),optional,intent(out) :: abs_n_tilde_nt_diff(nspden) 122 real(dp),optional,intent(in) :: znucl(ntypat) 123 type(pawrad_type),intent(in) :: pawrad(ntypat) 124 type(pawrhoij_type),intent(in) :: pawrhoij(my_natom) 125 type(pawtab_type),target,intent(in) :: pawtab(ntypat) 126 127 !Local variables------------------------------- 128 !scalars 129 integer,parameter :: master=0 130 integer :: delta,iatom,ierr,ifgd,ifftsph,inl,inrm,ipsang,irhoij 131 integer :: ispden,itypat,il,im,ilm,iln,ilmn 132 integer :: jl,jlm,jln,jm,j0lmn,jlmn 133 integer :: klmn,my_comm_atom,my_start_indx,my_end_indx 134 integer :: nfgd,nnl,normchoice,nprocs,optcut,optgr0,optgr1,optgr2 135 integer :: optrad,option,my_rank,remainder,tmp_unt 136 real(dp) :: phj,phi,rR,tphj,tphi,ybcbeg,ybcend 137 logical :: my_atmtab_allocated,paral_atom 138 character(len=500) :: message 139 !arrays 140 integer,allocatable :: l_size_atm(:),nrm_ifftsph(:) 141 integer,ABI_CONTIGUOUS pointer :: indlmn(:,:) 142 integer,pointer :: my_atmtab(:) 143 real(dp) :: ylmgr(3,3,0) 144 real(dp) :: yvals(4),xcart(3,natom) 145 real(dp),allocatable :: diag(:),nrm(:),phigrd(:,:),tphigrd(:,:),ylm(:,:),ypp(:) 146 real(dp),allocatable :: phi_at_zero(:),tphi_at_zero(:) 147 real(dp),allocatable :: rhor_tmp(:,:),tot_rhor(:) 148 character(len=fnlen) :: xsf_fname 149 type(pawfgrtab_type) :: local_pawfgrtab(my_natom) 150 151 ! ************************************************************************ 152 153 DBG_ENTER("COLL") 154 155 !Set up parallelism over atoms (compatible only with band-FFT parallelism) 156 paral_atom=(present(comm_atom).and.(my_natom/=natom)) 157 nullify(my_atmtab);if (present(mpi_atmtab)) my_atmtab => mpi_atmtab 158 my_comm_atom=xmpi_comm_self;if (present(comm_atom)) my_comm_atom=comm_atom 159 call get_my_atmtab(my_comm_atom,my_atmtab,my_atmtab_allocated,paral_atom,natom,& 160 & my_natom_ref=my_natom) 161 162 !MG: FIXME It won't work if atom-parallelism is used 163 !but even the loop over atoms below should be rewritten in this case. 164 165 !use a local copy of pawfgrtab to make sure we use the correction in the paw spheres 166 !the usual pawfgrtab uses r_shape which may not be the same as r_paw 167 if (my_natom>0) then 168 if (paral_atom) then 169 call pawtab_get_lsize(pawtab,l_size_atm,my_natom,typat,mpi_atmtab=my_atmtab) 170 call pawfgrtab_init(local_pawfgrtab,pawrhoij(1)%cplex,l_size_atm,nspden,typat,& 171 & mpi_atmtab=my_atmtab,comm_atom=my_comm_atom) 172 else 173 call pawtab_get_lsize(pawtab,l_size_atm,my_natom,typat) 174 call pawfgrtab_init(local_pawfgrtab,pawrhoij(1)%cplex,l_size_atm,nspden,typat) 175 end if 176 ABI_DEALLOCATE(l_size_atm) 177 end if 178 179 !Note: call to nhatgrid: comm_fft not used because FFT parallelism 180 !is done manually below 181 optcut = 1 ! use rpaw to construct local_pawfgrtab 182 optgr0 = 0; optgr1 = 0; optgr2 = 0 ! dont need gY terms locally 183 optrad = 1 ! do store r-R 184 if (paral_atom) then 185 call nhatgrid(atindx1,gmet,my_natom,natom,nattyp,ngfft,ntypat,& 186 & optcut,optgr0,optgr1,optgr2,optrad,local_pawfgrtab,pawtab,rprimd,typat,ucvol,xred,& 187 & comm_atom=my_comm_atom,mpi_atmtab=my_atmtab) 188 else 189 call nhatgrid(atindx1,gmet,my_natom,natom,nattyp,ngfft,ntypat,& 190 & optcut,optgr0,optgr1,optgr2,optrad,local_pawfgrtab,pawtab,rprimd,typat,ucvol,xred) 191 end if 192 !now local_pawfgrtab is ready to use 193 194 !Initialise output arrays. 195 rhor_paw=zero; rhor_n_one=zero; rhor_nt_one=zero 196 197 !Initialise and check parallell execution 198 my_rank = xmpi_comm_rank(spaceComm_in) 199 nprocs = xmpi_comm_size(spaceComm_in) 200 201 202 !loop over atoms in cell 203 do iatom = 1, my_natom 204 itypat = pawrhoij(iatom)%itypat 205 indlmn => pawtab(itypat)%indlmn 206 nfgd = local_pawfgrtab(iatom)%nfgd ! number of points in the fine grid for this PAW sphere 207 nnl = pawtab(itypat)%basis_size ! number of nl elements in PAW basis 208 209 ! Division of fine grid points among processors 210 if (nprocs==1) then ! Make sure everything runs with one proc 211 write(message,'(a)') ' In denfgr - number of processors: 1' 212 call wrtout(std_out,message,'COLL') 213 write(message,'(a)') ' Calculation of PAW density done in serial' 214 call wrtout(std_out,message,'COLL') 215 write(message,'(a,I6)') ' Number of fine grid points:',nfgd 216 call wrtout(std_out,message,'COLL') 217 my_start_indx = 1 218 my_end_indx = nfgd 219 else ! Divide up the fine grid points among the processors 220 write(message,'(a,I4)') ' In denfgr - number of processors: ',nprocs 221 call wrtout(std_out,message,'COLL') 222 write(message,'(a)') ' Calculation of PAW density done in parallel' 223 call wrtout(std_out,message,'COLL') 224 write(message,'(a,I6)') ' Number of fine grid points:',nfgd 225 call wrtout(std_out,message,'COLL') 226 ! Divide the fine grid points among the processors 227 delta = int(floor(real(nfgd)/real(nprocs))) 228 remainder = nfgd-nprocs*delta 229 my_start_indx = 1+my_rank*delta 230 my_end_indx = (my_rank+1)*delta 231 ! Divide the remainder points among the processors 232 ! by shuffling indices 233 if ((my_rank+1)>remainder) then 234 my_start_indx = my_start_indx + remainder 235 my_end_indx = my_end_indx + remainder 236 else 237 my_start_indx = my_start_indx + my_rank 238 my_end_indx = my_end_indx + my_rank + 1 239 end if 240 if (prtvol>9) then 241 write(message,'(a,I6)') ' My index Starts at: ',my_start_indx 242 call wrtout(std_out,message,'PERS') 243 write(message,'(a,I6)') ' Ends at: ',my_end_indx 244 call wrtout(std_out,message,'PERS') 245 write(message,'(a,I6)') ' # pts: ',my_end_indx+1-my_start_indx 246 call wrtout(std_out,message,'PERS') 247 end if 248 end if 249 250 write(message,'(a,I3,a,I3)') ' Entered loop for atom: ',iatom,' of:',natom 251 call wrtout(std_out,message,'PERS') 252 253 ! obtain |r-R| values on fine grid 254 ABI_ALLOCATE(nrm,(nfgd)) 255 do ifgd=1, nfgd 256 nrm(ifgd) = sqrt(dot_product(local_pawfgrtab(iatom)%rfgd(:,ifgd),local_pawfgrtab(iatom)%rfgd(:,ifgd))) 257 end do ! these are the |r-R| values 258 259 ! compute Ylm for each r-R vector. 260 ! ---- 261 ipsang = 1 + (pawtab(itypat)%l_size - 1)/2 ! recall l_size=2*l_max+1 262 ABI_ALLOCATE(ylm,(ipsang*ipsang,nfgd)) 263 option = 1 ! compute Ylm(r-R) for vectors 264 normchoice = 1 ! use computed norms of input vectors 265 call initylmr(ipsang,normchoice,nfgd,nrm,option,local_pawfgrtab(iatom)%rfgd,ylm,ylmgr) 266 267 ! in order to do spline fits, the |r-R| data must be sorted 268 ! ---- 269 ABI_ALLOCATE(nrm_ifftsph,(nfgd)) 270 nrm_ifftsph(:) = local_pawfgrtab(iatom)%ifftsph(:) ! copy of indices of points, to be rearranged by sort_dp 271 call sort_dp(nfgd,nrm,nrm_ifftsph,tol8) ! sort the nrm points, keeping track of which goes where 272 273 ! now make spline fits of phi and tphi onto the fine grid around the atom 274 ! ---- 275 ABI_ALLOCATE(phigrd,(nfgd,nnl)) 276 ABI_ALLOCATE(tphigrd,(nfgd,nnl)) 277 ABI_ALLOCATE(phi_at_zero,(nnl)) 278 ABI_ALLOCATE(tphi_at_zero,(nnl)) 279 ABI_ALLOCATE(ypp,(pawtab(itypat)%mesh_size)) 280 ABI_ALLOCATE(diag,(pawtab(itypat)%mesh_size)) 281 282 do inl = 1, nnl 283 284 ! spline phi onto points 285 ypp(:) = zero; diag(:) = zero; ybcbeg = zero; ybcend = zero; 286 call spline(pawrad(itypat)%rad,pawtab(itypat)%phi(:,inl),pawtab(itypat)%mesh_size,ybcbeg,ybcend,ypp) 287 call splint(pawtab(itypat)%mesh_size,pawrad(itypat)%rad,pawtab(itypat)%phi(:,inl),ypp,nfgd,nrm,phigrd(:,inl)) 288 289 ! next splint tphi onto points 290 ypp(:) = zero; diag(:) = zero; ybcbeg = zero; ybcend = zero; 291 call spline(pawrad(itypat)%rad,pawtab(itypat)%tphi(:,inl),pawtab(itypat)%mesh_size,ybcbeg,ybcend,ypp) 292 call splint(pawtab(itypat)%mesh_size,pawrad(itypat)%rad,pawtab(itypat)%tphi(:,inl),ypp,nfgd,nrm,tphigrd(:,inl)) 293 294 ! Find out the value of the basis function at zero using extrapolation 295 yvals = zero 296 ! Extrapolate only if this is an s-state (l=0) 297 if (indlmn(1,inl)==0) then 298 yvals(2:4) = pawtab(itypat)%phi(2:4,inl)/pawrad(itypat)%rad(2:4) 299 call pawrad_deducer0(yvals,4,pawrad(itypat)) 300 write(std_out,*) 'phi_at_zero: ',yvals(1),' from:',yvals(2:4) 301 end if 302 phi_at_zero(inl) = yvals(1) 303 304 yvals = zero 305 ! Extrapolate only if this is an s-state (l=0) 306 if (indlmn(1,inl)==0) then 307 yvals(2:4) = pawtab(itypat)%tphi(2:4,inl)/pawrad(itypat)%rad(2:4) 308 call pawrad_deducer0(yvals,4,pawrad(itypat)) 309 write(std_out,*) 'tphi_at_zero: ',yvals(1),' from:',yvals(2:4) 310 end if 311 tphi_at_zero(inl) = yvals(1) 312 313 end do ! end loop over nnl basis functions 314 ABI_DEALLOCATE(ypp) 315 ABI_DEALLOCATE(diag) 316 317 ! loop over basis elements for this atom 318 ! because we have to store things like <phi|r'><r'|phi>-<tphi|r'><r'|tphi> at each point of the 319 ! fine grid, there is no integration, and hence no simplifications of the Y_lm's. That's why 320 ! we have to loop through the basis elements in exhaustive detail, rather than just a loop over 321 ! lmn2_size or something comparable. 322 ! ---- 323 if (prtvol>9) then 324 write(message,'(a,I3)') ' Entering j-loop over basis elements for atom:',iatom 325 call wrtout(std_out,message,'PERS') 326 end if 327 328 do jlmn=1,pawtab(itypat)%lmn_size 329 330 if (prtvol>9) then 331 write(message,'(2(a,I3))') ' Element:',jlmn,' of:',pawtab(itypat)%lmn_size 332 call wrtout(std_out,message,'PERS') 333 end if 334 335 jl=indlmn(1,jlmn) 336 jm=indlmn(2,jlmn) 337 jlm=indlmn(4,jlmn) 338 jln=indlmn(5,jlmn) 339 j0lmn=jlmn*(jlmn-1)/2 340 341 if (prtvol>9) then 342 write(message,'(a,I3)') ' Entering i-loop for j:',jlmn 343 call wrtout(std_out,message,'PERS') 344 end if 345 346 do ilmn=1,jlmn 347 348 if (prtvol>9) then 349 write(message,'(2(a,I3))') ' Element:',ilmn,' of:',jlmn 350 call wrtout(std_out,message,'PERS') 351 end if 352 353 il=indlmn(1,ilmn) 354 im=indlmn(2,ilmn) 355 iln=indlmn(5,ilmn) 356 ilm=indlmn(4,ilmn) 357 klmn=j0lmn+ilmn 358 359 if (prtvol>9) then 360 write(message,'(a)') ' Entering loop over nonzero elems of rhoij' 361 call wrtout(std_out,message,'PERS') 362 end if 363 364 ! Loop over non-zero elements of rhoij 365 do irhoij=1,pawrhoij(iatom)%nrhoijsel 366 if (klmn==pawrhoij(iatom)%rhoijselect(irhoij)) then ! rho_ij /= 0 for this klmn 367 368 do ifgd=my_start_indx, my_end_indx ! loop over fine grid points in current PAW sphere 369 ifftsph = local_pawfgrtab(iatom)%ifftsph(ifgd) ! index of the point on the grid 370 371 ! have to retrieve the spline point to use since these were sorted 372 do inrm=1, nfgd 373 if(nrm_ifftsph(inrm) == ifftsph) exit ! have found nrm point corresponding to nfgd point 374 end do ! now inrm is the index of the sorted nrm vector to use 375 376 ! avoid division by zero 377 if(nrm(inrm) > zero) then 378 rR = nrm(inrm) ! value of |r-R| in the following 379 ! recall that <r|phi>=u(r)*Slm(r^)/r 380 phj = phigrd(inrm,jln)*ylm(jlm,ifgd)/rR 381 phi = phigrd(inrm,iln)*ylm(ilm,ifgd)/rR 382 tphj = tphigrd(inrm,jln)*ylm(jlm,ifgd)/rR 383 tphi = tphigrd(inrm,iln)*ylm(ilm,ifgd)/rR 384 else 385 ! use precalculated <r|phi>=u(r)*Slm(r^)/r at r=0 386 phj = phi_at_zero(jln)*ylm(jlm,ifgd) 387 phi = phi_at_zero(iln)*ylm(ilm,ifgd) 388 tphj = tphi_at_zero(jln)*ylm(jlm,ifgd) 389 tphi = tphi_at_zero(iln)*ylm(ilm,ifgd) 390 end if ! check if |r-R| = 0 391 392 do ispden=1,nspden 393 if (pawrhoij(iatom)%cplex == 1) then 394 rhor_paw(ifftsph,ispden) = rhor_paw(ifftsph,ispden) + & 395 & pawtab(itypat)%dltij(klmn)*pawrhoij(iatom)%rhoijp(irhoij,ispden)*(phj*phi - tphj*tphi) 396 397 rhor_n_one(ifftsph,ispden) = rhor_n_one(ifftsph,ispden) + & 398 & pawtab(itypat)%dltij(klmn)*pawrhoij(iatom)%rhoijp(irhoij,ispden)*phj*phi 399 400 rhor_nt_one(ifftsph,ispden) = rhor_nt_one(ifftsph,ispden) + & 401 & pawtab(itypat)%dltij(klmn)*pawrhoij(iatom)%rhoijp(irhoij,ispden)*tphj*tphi 402 else 403 rhor_paw(ifftsph,ispden) = rhor_paw(ifftsph,ispden) + & 404 & pawtab(itypat)%dltij(klmn)*pawrhoij(iatom)%rhoijp(2*irhoij-1,ispden)*(phj*phi - tphj*tphi) 405 406 rhor_n_one(ifftsph,ispden) = rhor_n_one(ifftsph,ispden) + & 407 & pawtab(itypat)%dltij(klmn)*pawrhoij(iatom)%rhoijp(2*irhoij-1,ispden)*phj*phi 408 409 rhor_nt_one(ifftsph,ispden) = rhor_nt_one(ifftsph,ispden) + & 410 & pawtab(itypat)%dltij(klmn)*pawrhoij(iatom)%rhoijp(2*irhoij-1,ispden)*tphj*tphi 411 end if ! end check on cplex rhoij 412 413 end do ! end loop over nsdpen 414 end do ! end loop over nfgd 415 end if ! end selection on rhoij /= 0 416 end do ! end loop over non-zero rhoij 417 end do ! end loop over ilmn atomic basis states 418 end do ! end loop over jlmn atomic basis states 419 420 ABI_DEALLOCATE(nrm) 421 ABI_DEALLOCATE(nrm_ifftsph) 422 ABI_DEALLOCATE(phigrd) 423 ABI_DEALLOCATE(tphigrd) 424 ABI_DEALLOCATE(ylm) 425 ABI_DEALLOCATE(phi_at_zero) 426 ABI_DEALLOCATE(tphi_at_zero) 427 end do ! Loop on atoms 428 429 !MPI sum on each node the different contributions to the PAW densities. 430 call xmpi_sum(rhor_paw,spaceComm_in,ierr) 431 call xmpi_sum(rhor_n_one,spaceComm_in,ierr) 432 call xmpi_sum(rhor_nt_one,spaceComm_in,ierr) 433 if (paral_atom) then 434 call xmpi_sum(rhor_paw,my_comm_atom,ierr) 435 call xmpi_sum(rhor_n_one,my_comm_atom,ierr) 436 call xmpi_sum(rhor_nt_one,my_comm_atom,ierr) 437 end if 438 439 call wrtout(std_out,' *** Partial contributions to PAW rhor summed ***','PERS') 440 call xmpi_barrier(spaceComm_in) 441 442 !Add the plane-wave contribution \tilde{n} and remove \hat{n} 443 !BE careful here since the storage mode of rhoij and rhor is different. 444 select case (nspinor) 445 case (1) 446 if (nsppol==1) then 447 rhor_paw = rhor_paw + rhor - nhat 448 else ! Spin-polarised case: rhor_paw contains rhor_paw(spin_up,spin_down) but we need rhor_paw(total,spin_up) 449 ABI_ALLOCATE(tot_rhor,(pawfgr%nfft)) 450 ! 451 ! AE rhor 452 tot_rhor(:) = SUM(rhor_paw,DIM=2) 453 rhor_paw(:,2) = rhor_paw(:,1) 454 rhor_paw(:,1) = tot_rhor 455 rhor_paw = rhor_paw + rhor - nhat 456 ! 457 ! onsite AE rhor 458 tot_rhor(:) = SUM(rhor_n_one,DIM=2) 459 rhor_n_one(:,2) = rhor_n_one(:,1) 460 rhor_n_one(:,1) = tot_rhor 461 ! 462 ! onsite PS rhor 463 tot_rhor(:) = SUM(rhor_nt_one,DIM=2) 464 rhor_nt_one(:,2) = rhor_nt_one(:,1) 465 rhor_nt_one(:,1) = tot_rhor 466 467 ABI_DEALLOCATE(tot_rhor) 468 end if 469 470 case (2) 471 ! * if nspden==4, rhor contains (n^11, n^22, Re[n^12], Im[n^12]. 472 ! Storage mode for rhoij is different, See pawaccrhoij. 473 MSG_ERROR("nspinor 2 not coded") 474 case default 475 write(message,'(a,i0)')" Wrong value for nspinor=",nspinor 476 MSG_ERROR(message) 477 end select 478 479 !if (prtvol>9) then ! Check normalisation 480 !write(message,'(a,F8.4)') ' PAWDEN - NORM OF DENSITY: ',SUM(rhor_paw(:,1))*ucvol/PRODUCT(pawfgr%ngfft(1:3)) 481 !call wrtout(std_out,message,'COLL') 482 !end if 483 484 if (present(abs_n_tilde_nt_diff).AND.present(znucl)) then 485 ABI_ALLOCATE(rhor_tmp,(pawfgr%nfft,nspden)) 486 do ispden=1,nspden 487 rhor_tmp(:,ispden) = zero 488 do iatom=1,my_natom 489 do ifgd=1,local_pawfgrtab(iatom)%nfgd ! loop over fine grid points in current PAW sphere 490 ifftsph = local_pawfgrtab(iatom)%ifftsph(ifgd) ! index of the point on the grid 491 rhor_tmp(ifftsph,ispden) = rhor(ifftsph,ispden) - nhat(ifftsph,ispden) & 492 & - rhor_nt_one(ifftsph,ispden) 493 end do !ifgd 494 end do ! iatom 495 end do ! ispden 496 if (paral_atom) then 497 call xmpi_sum(rhor_tmp,my_comm_atom,ierr) 498 end if 499 500 if (my_rank==master) then 501 do ispden=1,nspden 502 ! Write to xsf file 503 call xred2xcart(natom,rprimd,xcart,xred) 504 write(xsf_fname,'(a,I0,a)') 'N_tilde_onsite_diff_sp',ispden,'.xsf' 505 if (open_file(xsf_fname,message, unit=tmp_unt,status='unknown',form='formatted') /= 0) then 506 MSG_ERROR(message) 507 end if 508 call printxsf(ngfft(1),ngfft(2),ngfft(3),rhor_tmp(:,ispden),rprimd,& 509 & (/zero,zero,zero/),natom,ntypat,typat,xcart,znucl,tmp_unt,0) 510 close(tmp_unt) 511 abs_n_tilde_nt_diff(ispden) = SUM(ABS(rhor_tmp(:,ispden)))/pawfgr%nfft 512 write(message,'(4(a),F16.9,2(a,I0),a)') ch10,' Wrote xsf file with \tilde{n}-\tilde{n}^1.',ch10,& 513 & ' Value of norm |\tilde{n}-\tilde{n}^1|:',& 514 & abs_n_tilde_nt_diff(ispden),' spin: ',ispden,' of ',nspden,ch10 515 call wrtout(std_out,message,'COLL') 516 end do 517 end if 518 ABI_DEALLOCATE(rhor_tmp) 519 520 else if ((present(abs_n_tilde_nt_diff).AND.(.NOT.present(znucl))) & 521 & .OR.(.NOT.present(abs_n_tilde_nt_diff).AND.(present(znucl)))) then 522 write(message,'(a)') ' Both abs_n_tilde_nt_diff *and* znucl must be passed',ch10,& 523 & 'to denfgr for |\tilde{n}-\tilde{n}^1| norm evaluation.' 524 MSG_ERROR(message) 525 end if 526 527 call pawfgrtab_free(local_pawfgrtab) 528 529 !Destroy atom table used for parallelism 530 call free_my_atmtab(my_atmtab,my_atmtab_allocated) 531 532 call xmpi_barrier(spaceComm_in) 533 534 DBG_EXIT("COLL") 535 536 end subroutine denfgr