Line data Source code
1 : !!****m* ABINIT/m_paw_mkrho
2 : !! NAME
3 : !! m_paw_mkrho
4 : !!
5 : !! FUNCTION
6 : !! This module contains routines used to compute PAW density on the real space fine grid.
7 : !!
8 : !! COPYRIGHT
9 : !! Copyright (C) 2018-2026 ABINIT group (MT, JWZ)
10 : !! This file is distributed under the terms of the
11 : !! GNU General Public License, see ~abinit/COPYING
12 : !! or http://www.gnu.org/copyleft/gpl.txt .
13 : !!
14 : !! SOURCE
15 :
16 : #if defined HAVE_CONFIG_H
17 : #include "config.h"
18 : #endif
19 :
20 : #include "abi_common.h"
21 :
22 : MODULE m_paw_mkrho
23 :
24 : use defs_basis
25 : use m_abicore
26 : use m_errors
27 : use m_xmpi
28 : use m_extfpmd
29 :
30 : use defs_abitypes, only : MPI_type
31 : use m_time, only : timab
32 : use m_pawang, only : pawang_type
33 : use m_pawrad, only : pawrad_type,pawrad_deducer0
34 : use m_pawtab, only : pawtab_type,pawtab_get_lsize
35 : use m_paw_sphharm, only : initylmr
36 : use m_pawfgrtab, only : pawfgrtab_type,pawfgrtab_init,pawfgrtab_free
37 : use m_pawrhoij, only : pawrhoij_type,pawrhoij_copy,pawrhoij_free_unpacked, &
38 : & pawrhoij_nullify,pawrhoij_free,pawrhoij_symrhoij
39 : use m_pawfgr, only : pawfgr_type
40 : use m_paw_nhat, only : pawmknhat,nhatgrid
41 : use m_paral_atom, only : get_my_atmtab,free_my_atmtab
42 : use m_fourier_interpol, only : transgrid
43 :
44 : use m_sort, only : sort_dp
45 : use m_splines, only : spline,splint
46 : use m_io_tools, only : open_file
47 : use m_geometry, only : xred2xcart
48 : use m_pptools, only : printxsf
49 : use m_fft, only : fourdp
50 :
51 : implicit none
52 :
53 : private
54 :
55 : !public procedures.
56 : public :: pawmkrho ! Build PAW electronic density on fine grid, including compensation charge density
57 : public :: denfgr ! Build complete PAW electronic density on fine grid, including on-site contributions
58 :
59 : CONTAINS !========================================================================================
60 : !!***
61 :
62 : !----------------------------------------------------------------------
63 :
64 : !!****f* m_paw_mkrho/pawmkrho
65 : !! NAME
66 : !! pawmkrho
67 : !!
68 : !! FUNCTION
69 : !! PAW only:
70 : !! Build total pseudo (compensated) density (\tild_rho + \hat_rho)
71 : !! Build compensation charge density (\hat_rho)
72 : !! Build occupation matrix (packed storage)
73 : !!
74 : !! INPUTS
75 : !! compute_rhor_rhog: if 1: set the computation of rhor and rhog in addition to the compensating charge.
76 : !! if 0: compute only the compensating charge
77 : !! cplex: if 1, real space 1-order functions on FFT grid are REAL, if 2, COMPLEX
78 : !! 1 for GS calculations
79 : !! extfpmd <type(extfpmd_type)>=--optional--extended first-principles molecular dynamics type
80 : !! gprimd(3,3)=dimensional primitive translations for reciprocal space(bohr^-1).
81 : !! indsym(4,nsym,natom)=indirect indexing array for atom labels
82 : !! ipert=index of perturbation if pawrhoij is a pertubed rhoij
83 : !! no meaning for ground-state calculations (should be 0)
84 : !! idir=direction of atomic displacement (in case of atomic displ. perturb.)
85 : !! mpi_enreg=information about MPI parallelization
86 : !! my_natom=number of atoms treated by current processor
87 : !! natom=number of atoms in cell
88 : !! nspden=number of spin-density components
89 : !! nsym=number of symmetry elements in space group
90 : !! ntypat=number of types of atoms in unit cell.
91 : !! paral_kgb=option for (kpt,g vectors,bands) parallelism
92 : !! pawang <type(pawang_type)>=angular mesh discretization and related data
93 : !! pawang_sym <type(pawang_type)>=angular data used for symmetrization
94 : !! optional parameter only needed for RF calculations
95 : !! pawfgr <type(paw_fgr_type)>=fine rectangular grid parameters
96 : !! pawfgrtab(natom) <type(pawfgrtab_type)>=atomic data given on fine rectangular grid
97 : !! pawprtvol=control print volume and debugging output for PAW
98 : !! pawrhoij0(natom) <type(pawrhoij_type)>= GS paw rhoij occupancies and related data (used only if ipert>0)
99 : !! optional parameter only needed for RF calculations
100 : !! pawtab(ntypat) <type(pawtab_type)>=paw tabulated starting data
101 : !! qphon(3)=wavevector of the phonon (RF only)
102 : !! rhopsg(2,pawfgr%nfftc)= pseudo density given on the coarse grid in reciprocal space
103 : !! rhopsr(pawfgr%nfftc,nspden)= pseudo density given on the coarse grid in real space
104 : !! rprimd(3,3)=real space primitive translations.
105 : !! symafm(nsym)=(anti)ferromagnetic part of symmetry operations
106 : !! symrec(3,3,nsym)=symmetries of group in terms of operations on
107 : !! reciprocal space primitive translations
108 : !! typat(natom)=type for each atom
109 : !! ucvol=volume of the unit cell
110 : !! xred(3,natom)= reduced atomic coordinates
111 : !!
112 : !! OUTPUT
113 : !! compch_fft=compensation charge inside spheres integrated over fine fft grid
114 : !! pawnhat(pawfgr%nfft,nspden)=compensation charge density on fine rectangular grid (optional argument)
115 : !! rhog(2,pawfgr%nfft)= compensated pseudo density given on the fine grid in reciprocal space
116 : !! This output is optional
117 : !! rhor(pawfgr%nfft,nspden)= compensated pseudo density given on the fine grid in real space
118 : !!
119 : !! SIDE EFFECTS
120 : !! pawrhoij(my_natom)= PAW occupancies
121 : !! At input : values at previous step in packed storage (pawrhoij()%rhoijp)
122 : !! At output: values (symmetrized) in packed storage (pawrhoij()%rhoijp)
123 : !! pawrhoij_unsym(:)= unsymmetrized PAW occupancies
124 : !! At input : values (unsymmetrized) in unpacked storage (pawrhoij()%rhoij_)
125 : !! At output: values in unpacked storage (pawrhoij()%rhoij_) are destroyed
126 : !!
127 : !! NOTES
128 : !! pawrhoij and pawrhoij_unsym can be identical (refer to the same pawrhoij datastructure).
129 : !! They should be different only if pawrhoij is distributed over atomic sites
130 : !! (in that case pawrhoij_unsym should not be distributed over atomic sites).
131 : !!
132 : !! SOURCE
133 :
134 44008 : subroutine pawmkrho(compute_rhor_rhog,compch_fft,cplex,gprimd,idir,indsym,ipert,mpi_enreg,&
135 22004 : & my_natom,natom,nspden,nsym,ntypat,paral_kgb,pawang,pawfgr,pawfgrtab,pawprtvol,&
136 44008 : & pawrhoij,pawrhoij_unsym,&
137 44008 : & pawtab,qphon,rhopsg,rhopsr,rhor,rprimd,symafm,symrec,typat,ucvol,usewvl,xred,&
138 27597 : & pawang_sym,pawnhat,pawnhatgr,pawrhoij0,rhog,extfpmd) ! optional arguments
139 :
140 : !Arguments ------------------------------------
141 : !scalars
142 : integer,intent(in) :: compute_rhor_rhog,cplex,idir,ipert,my_natom,natom,nspden,nsym,ntypat,paral_kgb,pawprtvol
143 : integer,intent(in) :: usewvl
144 : real(dp),intent(in) :: ucvol
145 : real(dp),intent(out) :: compch_fft
146 : type(MPI_type),intent(in) :: mpi_enreg
147 : type(pawang_type),intent(in) :: pawang
148 : type(pawang_type),intent(in),optional :: pawang_sym
149 : type(pawfgr_type),intent(in) :: pawfgr
150 : type(extfpmd_type),intent(in),pointer,optional :: extfpmd
151 : !arrays
152 : integer,intent(in) :: indsym(4,nsym,natom)
153 : integer,intent(in) :: symafm(nsym),symrec(3,3,nsym),typat(natom)
154 : real(dp),intent(in) :: gprimd(3,3),qphon(3),rprimd(3,3),xred(3,natom)
155 : real(dp),intent(inout),target,optional :: pawnhat(cplex*pawfgr%nfft,nspden) !vz_i
156 : real(dp),intent(inout),target,optional :: pawnhatgr(:,:,:) !vz_i
157 : real(dp),intent(inout) :: rhor(cplex*pawfgr%nfft,nspden*compute_rhor_rhog)
158 : real(dp),intent(out),optional :: rhog(2,pawfgr%nfft*compute_rhor_rhog)
159 : real(dp),intent(inout) :: rhopsg(2,pawfgr%nfftc*compute_rhor_rhog),rhopsr(cplex*pawfgr%nfftc,nspden*compute_rhor_rhog)
160 : type(pawfgrtab_type),intent(inout) :: pawfgrtab(my_natom)
161 : type(pawrhoij_type),intent(inout),target :: pawrhoij(:)
162 : type(pawrhoij_type),intent(inout) :: pawrhoij_unsym(:)
163 : type(pawrhoij_type),intent(in),target,optional :: pawrhoij0(my_natom)
164 : type(pawtab_type),intent(in) :: pawtab(ntypat)
165 :
166 : !Local variables-------------------------------
167 : !scalars
168 : integer :: choice,ider,izero,option
169 : character(len=500) :: msg
170 : !arrays
171 : real(dp) :: tsec(2)
172 : real(dp),target :: rhodum(0,0,0)
173 22004 : real(dp),pointer :: pawnhat_ptr(:,:)
174 22004 : real(dp),pointer :: pawnhatgr_ptr(:,:,:)
175 22004 : type(pawrhoij_type),pointer :: pawrhoij_ptr(:),pawrhoij0_ptr(:)
176 :
177 : ! ***********************************************************************
178 :
179 : DBG_ENTER("COLL")
180 :
181 22004 : call timab(556,1,tsec)
182 :
183 : !Compatibility tests
184 22004 : if (size(pawrhoij_unsym)>0) then
185 22004 : if (pawrhoij_unsym(1)%use_rhoij_==0) then
186 0 : msg=' rhoij_ field must be allocated in pawrhoij_unsym !'
187 0 : ABI_BUG(msg)
188 : end if
189 : end if
190 22004 : if (ipert>0.and.(.not.present(pawrhoij0))) then
191 0 : msg=' pawrhoij0 must be present when ipert>0 !'
192 0 : ABI_BUG(msg)
193 : end if
194 :
195 : !Symetrize PAW occupation matrix and store it in packed storage
196 22004 : call timab(557,1,tsec)
197 22004 : option=1;choice=1
198 22004 : if (present(pawang_sym)) then
199 : call pawrhoij_symrhoij(pawrhoij,pawrhoij_unsym,choice,gprimd,indsym,ipert,&
200 : & natom,nsym,ntypat,option,pawang_sym,pawprtvol,pawtab,rprimd,symafm,&
201 : & symrec,typat,comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab,&
202 5708 : & qphon=qphon)
203 : else
204 : call pawrhoij_symrhoij(pawrhoij,pawrhoij_unsym,choice,gprimd,indsym,ipert,&
205 : & natom,nsym,ntypat,option,pawang,pawprtvol,pawtab,rprimd,symafm,&
206 : & symrec,typat,comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab,&
207 16296 : & qphon=qphon)
208 : end if
209 22004 : call pawrhoij_free_unpacked(pawrhoij_unsym)
210 22004 : call timab(557,2,tsec)
211 :
212 : !In somes cases (parallelism), has to distribute the PAW occupation matrix
213 22004 : if (size(pawrhoij)==natom.and.(my_natom/=natom)) then
214 0 : ABI_MALLOC(pawrhoij_ptr,(my_natom))
215 0 : call pawrhoij_nullify(pawrhoij_ptr)
216 : call pawrhoij_copy(pawrhoij,pawrhoij_ptr,&
217 : & mpi_atmtab=mpi_enreg%my_atmtab,comm_atom=mpi_enreg%comm_atom, &
218 0 : & keep_cplex=.false.,keep_qphase=.false.,keep_itypat=.false.,keep_nspden=.false.)
219 : else
220 22004 : pawrhoij_ptr=>pawrhoij
221 : end if
222 :
223 : !Compute compensation charge density
224 22004 : ider=0;izero=0
225 22004 : if (present(pawnhat)) then
226 21906 : pawnhat_ptr => pawnhat
227 : else
228 392 : ABI_MALLOC(pawnhat_ptr,(pawfgr%nfft,nspden))
229 : end if
230 22004 : if (present(pawnhatgr)) then
231 0 : pawnhatgr_ptr => pawnhatgr
232 : else
233 22004 : pawnhatgr_ptr => rhodum
234 : end if
235 22004 : if (present(pawrhoij0)) then
236 5593 : pawrhoij0_ptr => pawrhoij0
237 : else
238 16411 : pawrhoij0_ptr => pawrhoij_ptr
239 : end if
240 :
241 : call pawmknhat(compch_fft,cplex,ider,idir,ipert,izero,gprimd,my_natom,natom,&
242 : & pawfgr%nfft,pawfgr%ngfft,ider,nspden,ntypat,pawang,pawfgrtab,&
243 : & pawnhatgr_ptr,pawnhat_ptr,pawrhoij_ptr,pawrhoij0_ptr,pawtab,qphon,rprimd,ucvol,usewvl,xred,&
244 : & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab,&
245 : & comm_fft=mpi_enreg%comm_fft,paral_kgb=paral_kgb,me_g0=mpi_enreg%me_g0,&
246 22004 : & distribfft=mpi_enreg%distribfft,mpi_comm_wvl=mpi_enreg%comm_wvl)
247 :
248 22004 : if (compute_rhor_rhog/=0) then
249 : ! Transfer pseudo density from coarse grid to fine grid
250 21884 : if(usewvl==0) then
251 21884 : call transgrid(cplex,mpi_enreg,nspden,+1,1,0,paral_kgb,pawfgr,rhopsg,rhodum,rhopsr,rhor)
252 : end if
253 :
254 : ! Add pseudo density and compensation charge density (on fine grid)
255 898926947 : rhor(:,:)=rhor(:,:)+pawnhat_ptr(:,:)
256 :
257 : ! Add extfpmd electrons contributions to density on fine grid.
258 21884 : if(present(extfpmd)) then
259 15109 : if(associated(extfpmd)) then
260 90 : if(extfpmd%version==10.and.allocated(extfpmd%nelectarr)) then
261 0 : rhor(:,:)=rhor(:,:)+extfpmd%nelectarr(:,:)/ucvol/nspden
262 : else
263 993332 : rhor(:,:)=rhor(:,:)+extfpmd%nelect/ucvol/nspden
264 : end if
265 : end if
266 : end if
267 :
268 : ! Compute compensated pseudo density in reciprocal space
269 21884 : if (present(rhog)) then
270 18943 : call fourdp(cplex,rhog,rhor(:,1),-1,mpi_enreg,pawfgr%nfft,1,pawfgr%ngfft,0)
271 : end if
272 : end if
273 :
274 : !Free temporary memory spaces
275 22004 : if (.not.present(pawnhat)) then
276 98 : ABI_FREE(pawnhat_ptr)
277 : end if
278 22004 : if (size(pawrhoij)==natom.and.(my_natom/=natom)) then
279 0 : call pawrhoij_free(pawrhoij_ptr)
280 0 : ABI_FREE(pawrhoij_ptr)
281 : end if
282 : nullify(pawnhat_ptr)
283 : nullify(pawnhatgr_ptr)
284 22004 : nullify(pawrhoij_ptr)
285 :
286 22004 : call timab(556,2,tsec)
287 :
288 : DBG_EXIT("COLL")
289 :
290 66012 : end subroutine pawmkrho
291 : !!***
292 :
293 : !----------------------------------------------------------------------
294 :
295 : !!****f* m_paw_mkrho/denfgr
296 : !! NAME
297 : !! denfgr
298 : !!
299 : !! FUNCTION
300 : !! Construct complete electron density on fine grid, by removing nhat
301 : !! and adding PAW corrections
302 : !!
303 : !! INPUTS
304 : !! atindx1(natom)=index table for atoms, inverse of atindx (see gstate.f)
305 : !! gmet(3,3)=reciprocal space metric tensor in bohr**-2.
306 : !! mpi_atmtab(:)=--optional-- indexes of the atoms treated by current proc
307 : !! comm_atom=--optional-- MPI communicator over atoms
308 : !! my_natom=number of atoms treated by current processor
309 : !! natom= number of atoms in cell
310 : !! nattyp(ntypat)= # atoms of each type.
311 : !! ngfft(18)=contain all needed information about 3D FFT (see NOTES at beginning of scfcv)
312 : !! nhat(pawfgr%nfft,nspden)= compensation charge density used in PAW
313 : !! nspinor=Number of spinor components
314 : !! nsppol=Number of independent spin components.
315 : !! nspden= number of spin densities
316 : !! ntypat= number of types of atoms in the cell
317 : !! pawfgr <type(pawfgr_type)>= data about the fine grid
318 : !! pawrad(ntypat) <type(pawrad_type)>= radial mesh data for each type of atom
319 : !! pawrhoij(natom) <type(pawrhoij_type)>= rho_ij data for each atom
320 : !! pawtab(ntypat) <type(pawtab_type)>= PAW functions around each type of atom
321 : !! rhor(pawfgr%nfft,nspden)= input density ($\tilde{n}+\hat{n}$ in PAW case)
322 : !! rprimd(3,3)=dimensional primitive translations for real space (bohr)
323 : !! typat(natom)= list of atom types
324 : !! ucvol=unit cell volume (bohr**3)
325 : !! xred(3,natom)=reduced dimensionless atomic coordinates
326 : !!
327 : !! OUTPUT
328 : !! rhor_paw(pawfgr%nfft,nspden)= full electron density on the fine grid
329 : !!
330 : !! NOTES
331 : !! In PAW calculations, the valence density present in rhor includes the
332 : !! compensation charge density $\hat{n}$, and also doesn't include the on-site
333 : !! PAW contributions. For post-processing and proper visualization it is necessary
334 : !! to use the full electronic density, which is what this subroutine constructs.
335 : !! Specifically, it removes $\hat{n}$ from rhor, and also computes the on-site PAW
336 : !! terms. This is nothing other than the proper PAW treatment of the density
337 : !! operator $|\mathbf{r}\rangle\langle\mathbf{r}|$, and yields the formula
338 : !! $$\tilde{n}+\sum_{ij}\rho_ij\left[\varphi_i(\mathbf{r})\varphi_j(\mathbf{r})-
339 : !! \tilde{\varphi}_i(\mathbf{r})\tilde{\varphi}_j(\mathbf{r})\right]$$
340 : !! Notice that this formula is expressed on the fine grid, and requires
341 : !! interpolating the PAW radial functions onto this grid, as well as calling
342 : !! initylmr in order to get the angular functions on the grid points.
343 : !!
344 : !! SOURCE
345 :
346 4 : subroutine denfgr(atindx1,gmet,spaceComm_in,my_natom,natom,nattyp,ngfft,nhat,nspinor,nsppol,nspden,ntypat, &
347 4 : & pawfgr,pawrad,pawrhoij,pawtab,prtvol,rhor,rhor_paw,rhor_n_one,rhor_nt_one,rprimd,typat,ucvol,xred,&
348 2 : & abs_n_tilde_nt_diff,znucl,mpi_atmtab,comm_atom) ! Optional arguments
349 :
350 : !Arguments ------------------------------------
351 : !scalars
352 : integer,intent(in) :: my_natom,natom,nspden,ntypat,prtvol,nsppol,nspinor
353 : integer,optional,intent(in) :: comm_atom
354 : real(dp),intent(in) :: ucvol
355 : type(pawfgr_type),intent(in) :: pawfgr
356 : !arrays
357 : integer,intent(in) :: spaceComm_in
358 : integer,intent(in) :: atindx1(natom),nattyp(ntypat),ngfft(18),typat(natom)
359 : integer,optional,target,intent(in) :: mpi_atmtab(:)
360 : real(dp),intent(in) :: gmet(3,3),nhat(pawfgr%nfft,nspden)
361 : real(dp),intent(in) :: rhor(pawfgr%nfft,nspden),rprimd(3,3)
362 : real(dp),intent(inout) :: xred(3,natom)
363 : real(dp),intent(out) :: rhor_paw(pawfgr%nfft,nspden)
364 : real(dp),intent(out) :: rhor_n_one(pawfgr%nfft,nspden)
365 : real(dp),intent(out) :: rhor_nt_one(pawfgr%nfft,nspden)
366 : real(dp),optional,intent(out) :: abs_n_tilde_nt_diff(nspden)
367 : real(dp),optional,intent(in) :: znucl(ntypat)
368 : type(pawrad_type),intent(in) :: pawrad(ntypat)
369 : type(pawrhoij_type),intent(in) :: pawrhoij(my_natom)
370 : type(pawtab_type),target,intent(in) :: pawtab(ntypat)
371 :
372 : !Local variables-------------------------------
373 : !scalars
374 : integer,parameter :: master=0
375 : integer :: delta,iatom,ierr,ifgd,ifftsph,inl,inrm,ipsang,irhoij
376 : integer :: ispden,itypat,il,im,ilm,iln,ilmn
377 : integer :: jl,jlm,jln,jm,j0lmn,jlmn
378 : integer :: klmn,my_comm_atom,my_start_indx,my_end_indx
379 : integer :: nfgd,nnl,normchoice,nprocs,optcut,optgr0,optgr1,optgr2
380 : integer :: optrad,option,my_rank,remainder,tmp_unt
381 : real(dp) :: phj,phi,rR,tphj,tphi,ybcbeg,ybcend
382 : logical :: my_atmtab_allocated,paral_atom
383 : character(len=500) :: message
384 : !arrays
385 4 : integer,allocatable :: l_size_atm(:),nrm_ifftsph(:)
386 2 : integer,ABI_CONTIGUOUS pointer :: indlmn(:,:)
387 2 : integer,pointer :: my_atmtab(:)
388 : real(dp) :: ylmgr(3,3,0)
389 4 : real(dp) :: yvals(4),xcart(3,natom)
390 4 : real(dp),allocatable :: diag(:),nrm(:),phigrd(:,:),tphigrd(:,:),ylm(:,:),ypp(:)
391 2 : real(dp),allocatable :: phi_at_zero(:),tphi_at_zero(:)
392 2 : real(dp),allocatable :: rhor_tmp(:,:),tot_rhor(:)
393 : character(len=fnlen) :: xsf_fname
394 50 : type(pawfgrtab_type) :: local_pawfgrtab(my_natom)
395 :
396 : ! ************************************************************************
397 :
398 : DBG_ENTER("COLL")
399 :
400 2 : if (my_natom>0) then
401 2 : ABI_CHECK(pawrhoij(1)%qphase==1,'denfgr not supposed to be called with qphase/=1!')
402 : end if
403 :
404 : !Set up parallelism over atoms (compatible only with band-FFT parallelism)
405 2 : paral_atom=(present(comm_atom).and.(my_natom/=natom))
406 2 : nullify(my_atmtab);if (present(mpi_atmtab)) my_atmtab => mpi_atmtab
407 2 : my_comm_atom=xmpi_comm_self;if (present(comm_atom)) my_comm_atom=comm_atom
408 : call get_my_atmtab(my_comm_atom,my_atmtab,my_atmtab_allocated,paral_atom,natom,&
409 2 : & my_natom_ref=my_natom)
410 :
411 : !MG: FIXME It won't work if atom-parallelism is used
412 : !but even the loop over atoms below should be rewritten in this case.
413 :
414 : !use a local copy of pawfgrtab to make sure we use the correction in the paw spheres
415 : !the usual pawfgrtab uses r_shape which may not be the same as r_paw
416 2 : if (my_natom>0) then
417 2 : if (paral_atom) then
418 0 : call pawtab_get_lsize(pawtab,l_size_atm,my_natom,typat,mpi_atmtab=my_atmtab)
419 : call pawfgrtab_init(local_pawfgrtab,pawrhoij(1)%qphase,l_size_atm,nspden,typat,&
420 0 : & mpi_atmtab=my_atmtab,comm_atom=my_comm_atom)
421 : else
422 2 : call pawtab_get_lsize(pawtab,l_size_atm,my_natom,typat)
423 2 : call pawfgrtab_init(local_pawfgrtab,pawrhoij(1)%qphase,l_size_atm,nspden,typat)
424 : end if
425 2 : ABI_FREE(l_size_atm)
426 : end if
427 :
428 : !Note: call to nhatgrid: comm_fft not used because FFT parallelism
429 : !is done manually below
430 2 : optcut = 1 ! use rpaw to construct local_pawfgrtab
431 2 : optgr0 = 0; optgr1 = 0; optgr2 = 0 ! dont need gY terms locally
432 2 : optrad = 1 ! do store r-R
433 2 : if (paral_atom) then
434 : call nhatgrid(atindx1,gmet,my_natom,natom,nattyp,ngfft,ntypat,&
435 : & optcut,optgr0,optgr1,optgr2,optrad,local_pawfgrtab,pawtab,rprimd,typat,ucvol,xred,&
436 0 : & comm_atom=my_comm_atom,mpi_atmtab=my_atmtab)
437 : else
438 : call nhatgrid(atindx1,gmet,my_natom,natom,nattyp,ngfft,ntypat,&
439 2 : & optcut,optgr0,optgr1,optgr2,optrad,local_pawfgrtab,pawtab,rprimd,typat,ucvol,xred)
440 : end if
441 : !now local_pawfgrtab is ready to use
442 :
443 : !Initialise output arrays.
444 208907 : rhor_paw=zero; rhor_n_one=zero; rhor_nt_one=zero
445 :
446 : !Initialise and check parallell execution
447 2 : my_rank = xmpi_comm_rank(spaceComm_in)
448 2 : nprocs = xmpi_comm_size(spaceComm_in)
449 :
450 :
451 : !loop over atoms in cell
452 6 : do iatom = 1, my_natom
453 4 : itypat = pawrhoij(iatom)%itypat
454 4 : indlmn => pawtab(itypat)%indlmn
455 4 : nfgd = local_pawfgrtab(iatom)%nfgd ! number of points in the fine grid for this PAW sphere
456 4 : nnl = pawtab(itypat)%basis_size ! number of nl elements in PAW basis
457 :
458 : ! Division of fine grid points among processors
459 4 : if (nprocs==1) then ! Make sure everything runs with one proc
460 4 : write(message,'(a)') ' In denfgr - number of processors: 1'
461 4 : call wrtout(std_out,message,'COLL')
462 4 : write(message,'(a)') ' Calculation of PAW density done in serial'
463 4 : call wrtout(std_out,message,'COLL')
464 4 : write(message,'(a,I6)') ' Number of fine grid points:',nfgd
465 4 : call wrtout(std_out,message,'COLL')
466 4 : my_start_indx = 1
467 4 : my_end_indx = nfgd
468 : else ! Divide up the fine grid points among the processors
469 0 : write(message,'(a,I4)') ' In denfgr - number of processors: ',nprocs
470 0 : call wrtout(std_out,message,'COLL')
471 0 : write(message,'(a)') ' Calculation of PAW density done in parallel'
472 0 : call wrtout(std_out,message,'COLL')
473 0 : write(message,'(a,I6)') ' Number of fine grid points:',nfgd
474 0 : call wrtout(std_out,message,'COLL')
475 : ! Divide the fine grid points among the processors
476 0 : delta = int(floor(real(nfgd)/real(nprocs)))
477 0 : remainder = nfgd-nprocs*delta
478 0 : my_start_indx = 1+my_rank*delta
479 0 : my_end_indx = (my_rank+1)*delta
480 : ! Divide the remainder points among the processors
481 : ! by shuffling indices
482 0 : if ((my_rank+1)>remainder) then
483 0 : my_start_indx = my_start_indx + remainder
484 0 : my_end_indx = my_end_indx + remainder
485 : else
486 0 : my_start_indx = my_start_indx + my_rank
487 0 : my_end_indx = my_end_indx + my_rank + 1
488 : end if
489 0 : if (prtvol>9) then
490 0 : write(message,'(a,I6)') ' My index Starts at: ',my_start_indx
491 0 : call wrtout(std_out,message,'PERS')
492 0 : write(message,'(a,I6)') ' Ends at: ',my_end_indx
493 0 : call wrtout(std_out,message,'PERS')
494 0 : write(message,'(a,I6)') ' # pts: ',my_end_indx+1-my_start_indx
495 0 : call wrtout(std_out,message,'PERS')
496 : end if
497 : end if
498 :
499 4 : write(message,'(a,I3,a,I3)') ' Entered loop for atom: ',iatom,' of:',natom
500 4 : call wrtout(std_out,message,'PERS')
501 :
502 : ! obtain |r-R| values on fine grid
503 12 : ABI_MALLOC(nrm,(nfgd))
504 6852 : do ifgd=1, nfgd
505 27396 : nrm(ifgd) = sqrt(dot_product(local_pawfgrtab(iatom)%rfgd(:,ifgd),local_pawfgrtab(iatom)%rfgd(:,ifgd)))
506 : end do ! these are the |r-R| values
507 :
508 : ! compute Ylm for each r-R vector.
509 : ! ----
510 4 : ipsang = 1 + (pawtab(itypat)%l_size - 1)/2 ! recall l_size=2*l_max+1
511 16 : ABI_MALLOC(ylm,(ipsang*ipsang,nfgd))
512 4 : option = 1 ! compute Ylm(r-R) for vectors
513 4 : normchoice = 1 ! use computed norms of input vectors
514 4 : call initylmr(ipsang,normchoice,nfgd,nrm,option,local_pawfgrtab(iatom)%rfgd,ylm,ylmgr)
515 :
516 : ! in order to do spline fits, the |r-R| data must be sorted
517 : ! ----
518 12 : ABI_MALLOC(nrm_ifftsph,(nfgd))
519 6852 : nrm_ifftsph(:) = local_pawfgrtab(iatom)%ifftsph(:) ! copy of indices of points, to be rearranged by sort_dp
520 4 : call sort_dp(nfgd,nrm,nrm_ifftsph,tol8) ! sort the nrm points, keeping track of which goes where
521 :
522 : ! now make spline fits of phi and tphi onto the fine grid around the atom
523 : ! ----
524 16 : ABI_MALLOC(phigrd,(nfgd,nnl))
525 12 : ABI_MALLOC(tphigrd,(nfgd,nnl))
526 12 : ABI_MALLOC(phi_at_zero,(nnl))
527 8 : ABI_MALLOC(tphi_at_zero,(nnl))
528 12 : ABI_MALLOC(ypp,(pawtab(itypat)%mesh_size))
529 8 : ABI_MALLOC(diag,(pawtab(itypat)%mesh_size))
530 :
531 20 : do inl = 1, nnl
532 :
533 : ! spline phi onto points
534 19856 : ypp(:) = zero; diag(:) = zero; ybcbeg = zero; ybcend = zero;
535 16 : call spline(pawrad(itypat)%rad,pawtab(itypat)%phi(:,inl),pawtab(itypat)%mesh_size,ybcbeg,ybcend,ypp)
536 16 : call splint(pawtab(itypat)%mesh_size,pawrad(itypat)%rad,pawtab(itypat)%phi(:,inl),ypp,nfgd,nrm,phigrd(:,inl))
537 :
538 : ! next splint tphi onto points
539 19856 : ypp(:) = zero; diag(:) = zero; ybcbeg = zero; ybcend = zero;
540 16 : call spline(pawrad(itypat)%rad,pawtab(itypat)%tphi(:,inl),pawtab(itypat)%mesh_size,ybcbeg,ybcend,ypp)
541 16 : call splint(pawtab(itypat)%mesh_size,pawrad(itypat)%rad,pawtab(itypat)%tphi(:,inl),ypp,nfgd,nrm,tphigrd(:,inl))
542 :
543 : ! Find out the value of the basis function at zero using extrapolation
544 16 : yvals = zero
545 : ! Extrapolate only if this is an s-state (l=0)
546 16 : if (indlmn(1,inl)==0) then
547 32 : yvals(2:4) = pawtab(itypat)%phi(2:4,inl)/pawrad(itypat)%rad(2:4)
548 8 : call pawrad_deducer0(yvals,4,pawrad(itypat))
549 8 : write(std_out,*) 'phi_at_zero: ',yvals(1),' from:',yvals(2:4)
550 : end if
551 16 : phi_at_zero(inl) = yvals(1)
552 :
553 16 : yvals = zero
554 : ! Extrapolate only if this is an s-state (l=0)
555 16 : if (indlmn(1,inl)==0) then
556 32 : yvals(2:4) = pawtab(itypat)%tphi(2:4,inl)/pawrad(itypat)%rad(2:4)
557 8 : call pawrad_deducer0(yvals,4,pawrad(itypat))
558 8 : write(std_out,*) 'tphi_at_zero: ',yvals(1),' from:',yvals(2:4)
559 : end if
560 20 : tphi_at_zero(inl) = yvals(1)
561 :
562 : end do ! end loop over nnl basis functions
563 4 : ABI_FREE(ypp)
564 4 : ABI_FREE(diag)
565 :
566 : ! loop over basis elements for this atom
567 : ! because we have to store things like <phi|r'><r'|phi>-<tphi|r'><r'|tphi> at each point of the
568 : ! fine grid, there is no integration, and hence no simplifications of the Y_lm's. That's why
569 : ! we have to loop through the basis elements in exhaustive detail, rather than just a loop over
570 : ! lmn2_size or something comparable.
571 : ! ----
572 4 : if (prtvol>9) then
573 0 : write(message,'(a,I3)') ' Entering j-loop over basis elements for atom:',iatom
574 0 : call wrtout(std_out,message,'PERS')
575 : end if
576 :
577 36 : do jlmn=1,pawtab(itypat)%lmn_size
578 :
579 32 : if (prtvol>9) then
580 0 : write(message,'(2(a,I3))') ' Element:',jlmn,' of:',pawtab(itypat)%lmn_size
581 0 : call wrtout(std_out,message,'PERS')
582 : end if
583 :
584 32 : jl=indlmn(1,jlmn)
585 32 : jm=indlmn(2,jlmn)
586 32 : jlm=indlmn(4,jlmn)
587 32 : jln=indlmn(5,jlmn)
588 32 : j0lmn=jlmn*(jlmn-1)/2
589 :
590 32 : if (prtvol>9) then
591 0 : write(message,'(a,I3)') ' Entering i-loop for j:',jlmn
592 0 : call wrtout(std_out,message,'PERS')
593 : end if
594 :
595 180 : do ilmn=1,jlmn
596 :
597 144 : if (prtvol>9) then
598 0 : write(message,'(2(a,I3))') ' Element:',ilmn,' of:',jlmn
599 0 : call wrtout(std_out,message,'PERS')
600 : end if
601 :
602 144 : il=indlmn(1,ilmn)
603 144 : im=indlmn(2,ilmn)
604 144 : iln=indlmn(5,ilmn)
605 144 : ilm=indlmn(4,ilmn)
606 144 : klmn=j0lmn+ilmn
607 :
608 144 : if (prtvol>9) then
609 0 : write(message,'(a)') ' Entering loop over nonzero elems of rhoij'
610 0 : call wrtout(std_out,message,'PERS')
611 : end if
612 :
613 : ! Loop over non-zero elements of rhoij
614 1904 : do irhoij=1,pawrhoij(iatom)%nrhoijsel
615 1872 : if (klmn==pawrhoij(iatom)%rhoijselect(irhoij)) then ! rho_ij /= 0 for this klmn
616 :
617 82224 : do ifgd=my_start_indx, my_end_indx ! loop over fine grid points in current PAW sphere
618 82176 : ifftsph = local_pawfgrtab(iatom)%ifftsph(ifgd) ! index of the point on the grid
619 :
620 : ! have to retrieve the spline point to use since these were sorted
621 113671320 : do inrm=1, nfgd
622 113671320 : if(nrm_ifftsph(inrm) == ifftsph) exit ! have found nrm point corresponding to nfgd point
623 : end do ! now inrm is the index of the sorted nrm vector to use
624 :
625 : ! avoid division by zero
626 82176 : if(nrm(inrm) > zero) then
627 82140 : rR = nrm(inrm) ! value of |r-R| in the following
628 : ! recall that <r|phi>=u(r)*Slm(r^)/r
629 82140 : phj = phigrd(inrm,jln)*ylm(jlm,ifgd)/rR
630 82140 : phi = phigrd(inrm,iln)*ylm(ilm,ifgd)/rR
631 82140 : tphj = tphigrd(inrm,jln)*ylm(jlm,ifgd)/rR
632 82140 : tphi = tphigrd(inrm,iln)*ylm(ilm,ifgd)/rR
633 : else
634 : ! use precalculated <r|phi>=u(r)*Slm(r^)/r at r=0
635 36 : phj = phi_at_zero(jln)*ylm(jlm,ifgd)
636 36 : phi = phi_at_zero(iln)*ylm(ilm,ifgd)
637 36 : tphj = tphi_at_zero(jln)*ylm(jlm,ifgd)
638 36 : tphi = tphi_at_zero(iln)*ylm(ilm,ifgd)
639 : end if ! check if |r-R| = 0
640 :
641 237720 : do ispden=1,nspden
642 237672 : if (pawrhoij(iatom)%cplex_rhoij == 1) then
643 : rhor_paw(ifftsph,ispden) = rhor_paw(ifftsph,ispden) + &
644 155496 : & pawtab(itypat)%dltij(klmn)*pawrhoij(iatom)%rhoijp(irhoij,ispden)*(phj*phi - tphj*tphi)
645 :
646 : rhor_n_one(ifftsph,ispden) = rhor_n_one(ifftsph,ispden) + &
647 155496 : & pawtab(itypat)%dltij(klmn)*pawrhoij(iatom)%rhoijp(irhoij,ispden)*phj*phi
648 :
649 : rhor_nt_one(ifftsph,ispden) = rhor_nt_one(ifftsph,ispden) + &
650 155496 : & pawtab(itypat)%dltij(klmn)*pawrhoij(iatom)%rhoijp(irhoij,ispden)*tphj*tphi
651 : else
652 : rhor_paw(ifftsph,ispden) = rhor_paw(ifftsph,ispden) + &
653 0 : & pawtab(itypat)%dltij(klmn)*pawrhoij(iatom)%rhoijp(2*irhoij-1,ispden)*(phj*phi - tphj*tphi)
654 :
655 : rhor_n_one(ifftsph,ispden) = rhor_n_one(ifftsph,ispden) + &
656 0 : & pawtab(itypat)%dltij(klmn)*pawrhoij(iatom)%rhoijp(2*irhoij-1,ispden)*phj*phi
657 :
658 : rhor_nt_one(ifftsph,ispden) = rhor_nt_one(ifftsph,ispden) + &
659 0 : & pawtab(itypat)%dltij(klmn)*pawrhoij(iatom)%rhoijp(2*irhoij-1,ispden)*tphj*tphi
660 : end if ! end check on cplex rhoij
661 :
662 : end do ! end loop over nsdpen
663 : end do ! end loop over nfgd
664 : end if ! end selection on rhoij /= 0
665 : end do ! end loop over non-zero rhoij
666 : end do ! end loop over ilmn atomic basis states
667 : end do ! end loop over jlmn atomic basis states
668 :
669 4 : ABI_FREE(nrm)
670 4 : ABI_FREE(nrm_ifftsph)
671 4 : ABI_FREE(phigrd)
672 4 : ABI_FREE(tphigrd)
673 4 : ABI_FREE(ylm)
674 4 : ABI_FREE(phi_at_zero)
675 6 : ABI_FREE(tphi_at_zero)
676 : end do ! Loop on atoms
677 :
678 : !MPI sum on each node the different contributions to the PAW densities.
679 2 : call xmpi_sum(rhor_paw,spaceComm_in,ierr)
680 2 : call xmpi_sum(rhor_n_one,spaceComm_in,ierr)
681 2 : call xmpi_sum(rhor_nt_one,spaceComm_in,ierr)
682 2 : if (paral_atom) then
683 0 : call xmpi_sum(rhor_paw,my_comm_atom,ierr)
684 0 : call xmpi_sum(rhor_n_one,my_comm_atom,ierr)
685 0 : call xmpi_sum(rhor_nt_one,my_comm_atom,ierr)
686 : end if
687 :
688 2 : call wrtout(std_out,' *** Partial contributions to PAW rhor summed ***','PERS')
689 2 : call xmpi_barrier(spaceComm_in)
690 :
691 : !Add the plane-wave contribution \tilde{n} and remove \hat{n}
692 : !BE careful here since the storage mode of rhoij and rhor is different.
693 4 : select case (nspinor)
694 : case (1)
695 2 : if (nsppol==1) then
696 4098 : rhor_paw = rhor_paw + rhor - nhat
697 : else ! Spin-polarised case: rhor_paw contains rhor_paw(spin_up,spin_down) but we need rhor_paw(total,spin_up)
698 3 : ABI_MALLOC(tot_rhor,(pawfgr%nfft))
699 : !
700 : ! AE rhor
701 98305 : tot_rhor(:) = SUM(rhor_paw,DIM=2)
702 32769 : rhor_paw(:,2) = rhor_paw(:,1)
703 32769 : rhor_paw(:,1) = tot_rhor
704 65539 : rhor_paw = rhor_paw + rhor - nhat
705 : !
706 : ! onsite AE rhor
707 98305 : tot_rhor(:) = SUM(rhor_n_one,DIM=2)
708 32769 : rhor_n_one(:,2) = rhor_n_one(:,1)
709 32769 : rhor_n_one(:,1) = tot_rhor
710 : !
711 : ! onsite PS rhor
712 98305 : tot_rhor(:) = SUM(rhor_nt_one,DIM=2)
713 32769 : rhor_nt_one(:,2) = rhor_nt_one(:,1)
714 32769 : rhor_nt_one(:,1) = tot_rhor
715 :
716 1 : ABI_FREE(tot_rhor)
717 : end if
718 :
719 : case (2)
720 : ! * if nspden==4, rhor contains (n^11, n^22, Re[n^12], Im[n^12].
721 : ! Storage mode for rhoij is different, See pawaccrhoij.
722 0 : ABI_ERROR("nspinor 2 not coded")
723 : case default
724 0 : write(message,'(a,i0)')" Wrong value for nspinor=",nspinor
725 2 : ABI_ERROR(message)
726 : end select
727 :
728 : !if (prtvol>9) then ! Check normalisation
729 : !write(message,'(a,F8.4)') ' PAWDEN - NORM OF DENSITY: ',SUM(rhor_paw(:,1))*ucvol/PRODUCT(pawfgr%ngfft(1:3))
730 : !call wrtout(std_out,message,'COLL')
731 : !end if
732 :
733 2 : if (present(abs_n_tilde_nt_diff).AND.present(znucl)) then
734 0 : ABI_MALLOC(rhor_tmp,(pawfgr%nfft,nspden))
735 0 : do ispden=1,nspden
736 0 : rhor_tmp(:,ispden) = zero
737 0 : do iatom=1,my_natom
738 0 : do ifgd=1,local_pawfgrtab(iatom)%nfgd ! loop over fine grid points in current PAW sphere
739 0 : ifftsph = local_pawfgrtab(iatom)%ifftsph(ifgd) ! index of the point on the grid
740 : rhor_tmp(ifftsph,ispden) = rhor(ifftsph,ispden) - nhat(ifftsph,ispden) &
741 0 : & - rhor_nt_one(ifftsph,ispden)
742 : end do !ifgd
743 : end do ! iatom
744 : end do ! ispden
745 0 : if (paral_atom) then
746 0 : call xmpi_sum(rhor_tmp,my_comm_atom,ierr)
747 : end if
748 :
749 0 : if (my_rank==master) then
750 0 : do ispden=1,nspden
751 : ! Write to xsf file
752 0 : call xred2xcart(natom,rprimd,xcart,xred)
753 0 : write(xsf_fname,'(a,I0,a)') 'N_tilde_onsite_diff_sp',ispden,'.xsf'
754 0 : if (open_file(xsf_fname,message, unit=tmp_unt,status='unknown',form='formatted') /= 0) then
755 0 : ABI_ERROR(message)
756 : end if
757 : call printxsf(ngfft(1),ngfft(2),ngfft(3),rhor_tmp(:,ispden),rprimd,&
758 0 : & (/zero,zero,zero/),natom,ntypat,typat,xcart,znucl,tmp_unt,0)
759 0 : close(tmp_unt)
760 0 : abs_n_tilde_nt_diff(ispden) = SUM(ABS(rhor_tmp(:,ispden)))/pawfgr%nfft
761 0 : write(message,'(4(a),F16.9,2(a,I0),a)') ch10,' Wrote xsf file with \tilde{n}-\tilde{n}^1.',ch10,&
762 0 : & ' Value of norm |\tilde{n}-\tilde{n}^1|:',&
763 0 : & abs_n_tilde_nt_diff(ispden),' spin: ',ispden,' of ',nspden,ch10
764 0 : call wrtout(std_out,message,'COLL')
765 : end do
766 : end if
767 0 : ABI_FREE(rhor_tmp)
768 :
769 : else if ((present(abs_n_tilde_nt_diff).AND.(.NOT.present(znucl))) &
770 2 : & .OR.(.NOT.present(abs_n_tilde_nt_diff).AND.(present(znucl)))) then
771 0 : write(message,'(a)') ' Both abs_n_tilde_nt_diff *and* znucl must be passed',ch10,&
772 0 : & 'to denfgr for |\tilde{n}-\tilde{n}^1| norm evaluation.'
773 0 : ABI_ERROR(message)
774 : end if
775 :
776 2 : call pawfgrtab_free(local_pawfgrtab)
777 :
778 : !Destroy atom table used for parallelism
779 2 : call free_my_atmtab(my_atmtab,my_atmtab_allocated)
780 :
781 2 : call xmpi_barrier(spaceComm_in)
782 :
783 : DBG_EXIT("COLL")
784 :
785 8 : end subroutine denfgr
786 : !!***
787 :
788 : !----------------------------------------------------------------------
789 :
790 : END MODULE m_paw_mkrho
791 : !!***
|
