Branch data 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 : 45776 : subroutine pawmkrho(compute_rhor_rhog,compch_fft,cplex,gprimd,idir,indsym,ipert,mpi_enreg,&
135 : 22888 : & my_natom,natom,nspden,nsym,ntypat,paral_kgb,pawang,pawfgr,pawfgrtab,pawprtvol,&
136 [ - + - + ]: 45776 : & pawrhoij,pawrhoij_unsym,&
137 : 45776 : & pawtab,qphon,rhopsg,rhopsr,rhor,rprimd,symafm,symrec,typat,ucvol,usewvl,xred,&
138 [ - + - - : 28545 : & 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 : 22888 : real(dp),pointer :: pawnhat_ptr(:,:)
174 : 22888 : real(dp),pointer :: pawnhatgr_ptr(:,:,:)
175 : 22888 : type(pawrhoij_type),pointer :: pawrhoij_ptr(:),pawrhoij0_ptr(:)
176 : :
177 : : ! ***********************************************************************
178 : :
179 : : DBG_ENTER("COLL")
180 : :
181 : 22888 : call timab(556,1,tsec)
182 : :
183 : : !Compatibility tests
184 [ + - ]: 22888 : if (size(pawrhoij_unsym)>0) then
185 [ - + ]: 22888 : 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 [ + + - + ]: 22888 : 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 : 22888 : call timab(557,1,tsec)
197 : 22888 : option=1;choice=1
198 [ + + ]: 22888 : 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 : 5772 : & 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 : 17116 : & qphon=qphon)
208 : : end if
209 : 22888 : call pawrhoij_free_unpacked(pawrhoij_unsym)
210 : 22888 : call timab(557,2,tsec)
211 : :
212 : : !In somes cases (parallelism), has to distribute the PAW occupation matrix
213 [ + + - + ]: 22888 : 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 : 22888 : pawrhoij_ptr=>pawrhoij
221 : : end if
222 : :
223 : : !Compute compensation charge density
224 : 22888 : ider=0;izero=0
225 [ + + ]: 22888 : if (present(pawnhat)) then
226 : 22790 : pawnhat_ptr => pawnhat
227 : : else
228 [ - + + - : 392 : ABI_MALLOC(pawnhat_ptr,(pawfgr%nfft,nspden))
- + + - +
- + - - +
- + ]
229 : : end if
230 [ - + ]: 22888 : if (present(pawnhatgr)) then
231 : 0 : pawnhatgr_ptr => pawnhatgr
232 : : else
233 : 22888 : pawnhatgr_ptr => rhodum
234 : : end if
235 [ + + ]: 22888 : if (present(pawrhoij0)) then
236 : 5657 : pawrhoij0_ptr => pawrhoij0
237 : : else
238 : 17231 : 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 [ - + - + : 22888 : & distribfft=mpi_enreg%distribfft,mpi_comm_wvl=mpi_enreg%comm_wvl)
- + - + ]
247 : :
248 [ + + ]: 22888 : if (compute_rhor_rhog/=0) then
249 : : ! Transfer pseudo density from coarse grid to fine grid
250 [ + - ]: 22768 : if(usewvl==0) then
251 : 22768 : 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 [ + + + + ]: 983516089 : rhor(:,:)=rhor(:,:)+pawnhat_ptr(:,:)
256 : :
257 : : ! Add extfpmd electrons contributions to density on fine grid.
258 [ + + ]: 22768 : if(present(extfpmd)) then
259 [ + + ]: 15541 : 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 [ + + ]: 22768 : if (present(rhog)) then
270 : 19803 : 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 [ + + ]: 22888 : if (.not.present(pawnhat)) then
276 [ - + ]: 98 : ABI_FREE(pawnhat_ptr)
277 : : end if
278 [ + + - + ]: 22888 : 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 : 22888 : nullify(pawrhoij_ptr)
285 : :
286 : 22888 : call timab(556,2,tsec)
287 : :
288 : : DBG_EXIT("COLL")
289 : :
290 [ + + + + ]: 68664 : 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,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 : : !!***
|
