Line data Source code
1 : !!****m* m_mkrho/m_mkrho
2 : !! NAME
3 : !! m_mkrho
4 : !!
5 : !! FUNCTION
6 : !! Procedures for computing densities from KS orbitals.
7 : !!
8 : !! COPYRIGHT
9 : !! Copyright (C) 1998-2026 ABINIT group (DCA, XG, GMR, LSI, AR, MB, MT, SM, VR, FJ)
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_mkrho
23 :
24 : use, intrinsic :: iso_c_binding, only: c_size_t, c_loc
25 :
26 : use defs_basis
27 : use defs_wvltypes
28 : use m_abicore
29 : use m_xmpi
30 : use m_xomp
31 : use m_errors
32 : use m_dtset
33 : use m_extfpmd
34 : use m_gputk
35 : use m_abi_linalg
36 :
37 : use defs_abitypes, only : MPI_type
38 : use m_fstrings, only : sjoin, itoa
39 : use m_time, only : timab
40 : use m_fftcore, only : sphereboundary
41 : use m_fft, only : fftpac, zerosym, fourwf, fourwf_optmem, fourdp
42 : use m_bandfft_kpt, only : bandfft_kpt_set_ikpt
43 : use m_paw_dmft, only : paw_dmft_type
44 : use m_spacepar, only : symrhg
45 : use defs_datatypes, only : pseudopotential_type
46 : use m_atomdata, only : atom_length
47 : use m_mpinfo, only : ptabs_fourdp, proc_distrb_cycle
48 : use m_pawtab, only : pawtab_type
49 : use m_io_tools, only : open_file
50 : use m_splines, only : spline, splint
51 : use m_sort, only : sort_dp
52 : use m_prep_kgb, only : prep_fourwf
53 : use m_wvl_rho, only : wvl_mkrho
54 : use m_rot_cg, only : rot_cg
55 :
56 : #if defined HAVE_YAKL
57 : use gator_mod
58 : #endif
59 :
60 : implicit none
61 :
62 : private
63 : !!***
64 :
65 : public :: mkrho
66 : public :: initro
67 : public :: prtrhomxmn
68 : public :: read_atomden
69 : !!***
70 :
71 : contains
72 : !!***
73 :
74 : !!****f* m_mkrho/mkrho
75 : !! NAME
76 : !! mkrho
77 : !!
78 : !! FUNCTION
79 : !! Depending on option argument value:
80 : !! --Compute charge density rho(r) and rho(G) in electrons/bohr**3
81 : !! from input wavefunctions, band occupations, and k point wts.
82 : !! --Compute kinetic energy density tau(r) and tau(G) in bohr**-5
83 : !! from input wavefunctions, band occupations, and k point wts.
84 : !! --Compute a given element of the kinetic energy density tensor
85 : !! tau_{alpha,beta}(r) and tau_{alpha,beta}(G) in bohr**-5
86 : !! from input wavefunctions, band occupations, and k point wts.
87 : !!
88 : !! INPUTS
89 : !! cg(2,mcg)=wf in G space
90 : !! dtset <type(dataset_type)>=all input variables for this dataset
91 : !! | istwfk(nkpt)=input option parameter that describes the storage of wfs
92 : !! | mband=maximum number of bands
93 : !! | mgfft=maximum size of 1D FFTs
94 : !! | mkmem=Number of k points treated by this node
95 : !! | mpw=maximum allowed value for npw
96 : !! | nband(nkpt*nsppol)=number of bands to be included in summation
97 : !! | at each k point for each spin channel
98 : !! | nfft=(effective) number of FFT grid points (for this processor)
99 : !! | ngfft(18)=contain all needed information about 3D FFT,
100 : !! | see ~abinit/doc/variables/vargs.htm#ngfft
101 : !! | nkpt=number of k points
102 : !! | nspden=number of spin-density components
103 : !! | nsppol=1 for unpolarized, 2 for spin-polarized
104 : !! | nsym=number of symmetry elements in group (at least 1 for identity)
105 : !! | symafm(nsym)=(anti)ferromagnetic part of symmetry operations
106 : !! | symrel(3,3,nsym)=symmetry matrices in real space (integers)
107 : !! | wtk(nkpt)=k point weights (they sum to 1.0)
108 : !! extfpmd <type(extfpmd_type)>=--optional--extended first-principles molecular dynamics type
109 : !! gprimd(3,3)=dimensional reciprocal space primitive translations
110 : !! irrzon(nfft**(1-1/nsym),2,(nspden/nsppol)-3*(nspden/4))=irreducible zone data
111 : !! kg(3,mpw*mkmem)=reduced planewave coordinates
112 : !! mcg=size of wave-functions array (cg) =mpw*nspinor*mband*mkmem*nsppol
113 : !! mpi_enreg=information about MPI parallelization
114 : !! npwarr(nkpt)=number of planewaves and boundary planewaves at each k
115 : !! occ(mband*nkpt*nsppol)=
116 : !! occupation numbers for each band (usually 2.0) at each k point
117 : !! option if 0: compute rhor (electron density)
118 : !! if 1: compute taur (kinetic energy density)
119 : !! (i.e. Trace over the kinetic energy density tensor)
120 : !! if 2: compute taur_{alpha,beta} !!NOT YET IMPLEMENTED
121 : !! (a given element of the kinetic energy density tensor)
122 : !! paw_dmft <type(paw_dmft_type)>= paw+dmft related data
123 : !! phnons(2,nfft**(1-1/nsym),(nspden/nsppol)-3*(nspden/4))=nonsymmorphic translation phases
124 : !! rprimd(3,3)=dimensional real space primitive translations
125 : !! tim_mkrho=timing code of the calling routine(can be set to 0 if not attributed)
126 : !! ucvol=unit cell volume (Bohr**3)
127 : !! wvl_den <type(wvl_denspot_type)>=density information for wavelets
128 : !! wvl_wfs <type(wvl_projector_type)>=wavefunctions information for wavelets
129 : !!
130 : !! OUTPUT
131 : !! rhog(2,nfft)=total electron density in G space
132 : !! rhor(nfft,nspden)=electron density in r space
133 : !! (if spin polarized, array contains total density in first half and spin-up density in second half)
134 : !! (for non-collinear magnetism, first element: total density, 3 next ones: mx,my,mz in units of hbar/2)
135 : !!
136 : !! SOURCE
137 :
138 17150 : subroutine mkrho(cg,dtset,gprimd,irrzon,kg,mcg,mpi_enreg,npwarr,occ,paw_dmft,phnons,&
139 17150 : & rhog,rhor,rprimd,tim_mkrho,ucvol,wvl_den,wvl_wfs,&
140 : & option,extfpmd,nfft_blocks, printout) !optional
141 :
142 : !Arguments ------------------------------------
143 : !scalars
144 : integer,intent(in) :: mcg,tim_mkrho
145 : integer,intent(in),optional :: option,nfft_blocks
146 : logical, optional :: printout ! Controls wether results are printed at the end. Default=yes.
147 : real(dp),intent(in) :: ucvol
148 : type(extfpmd_type),intent(in),pointer,optional :: extfpmd
149 : type(MPI_type),intent(in) :: mpi_enreg
150 : type(dataset_type),intent(in) :: dtset
151 : type(paw_dmft_type), intent(in) :: paw_dmft
152 : type(wvl_wf_type),intent(inout) :: wvl_wfs
153 : type(wvl_denspot_type), intent(inout) :: wvl_den
154 : !nfft**(1-1/nsym) is 1 if nsym==1, and nfft otherwise
155 : integer, intent(in) :: irrzon(dtset%nfft**(1-1/dtset%nsym),2, (dtset%nspden/dtset%nsppol)-3*(dtset%nspden/4))
156 : integer, intent(in) :: kg(3,dtset%mpw*dtset%mkmem),npwarr(dtset%nkpt)
157 : real(dp), intent(in) :: gprimd(3,3)
158 : real(dp), intent(in), target :: cg(2,mcg)
159 : real(dp), intent(in) :: occ(dtset%mband*dtset%nkpt*dtset%nsppol)
160 : !nfft**(1-1/nsym) is 1 if nsym==1, and nfft otherwise
161 : real(dp), intent(in) :: phnons(2,(dtset%ngfft(1)*dtset%ngfft(2)*dtset%ngfft(3))**(1-1/dtset%nsym), &
162 : (dtset%nspden/dtset%nsppol)-3*(dtset%nspden/4))
163 : real(dp), intent(in) :: rprimd(3,3)
164 : real(dp), intent(out) :: rhor(dtset%nfft,dtset%nspden),rhog(2,dtset%nfft)
165 :
166 : !Local variables-------------------------------
167 : !scalars
168 : integer,save :: nskip=0
169 : integer :: alpha,use_nondiag_occup_dmft,bdtot_index,beta,blocksize,iband,iband1,ibandc1,ib,iblock,icg,ierr
170 : integer :: ifft,ikg,ikpt,ioption,ipw,ipwbd,ipwsp,ishf,ispden,ispinor,ispinor_index
171 : integer :: isppol,istwf_k,jspinor_index
172 : integer :: me,my_nspinor,n1,n2,n3,n4,n5,n6,nalpha,nband_k,nband_occ,nbandc1,nbdblock,nbeta
173 : integer :: ndat,nfftot,npw_k,spaceComm,tim_fourwf,gpu_option,l_nfft_blocks,nfft_blocks_occ,nband_fftblock_occ
174 : integer :: iband_me
175 : integer :: mband_mem
176 : logical :: transfer_cg
177 : real(dp) :: kpt_cart,kg_k_cart,gp2pi1,gp2pi2,gp2pi3,cwftmp
178 : real(dp) :: weight,weight_i
179 : logical :: l_printout
180 : !character(len=500) :: message
181 : !arrays
182 17150 : integer,allocatable :: gbound(:,:)
183 : integer, contiguous, pointer :: kg_k(:,:) => null()
184 : logical :: locc_test,nspinor1TreatedByThisProc,nspinor2TreatedByThisProc,gpu_cwavef
185 : real(dp) :: dummy(2,1) = reshape( (/0.0, 0.0/), shape(dummy))
186 : real(dp) :: tsec(2)
187 17150 : real(dp),allocatable :: cwavef_rot(:,:,:,:),occ_diag(:),occ_k(:)
188 17150 : real(dp),allocatable :: kg_k_cart_block(:),taur_alphabeta(:,:,:,:),weight_t(:)
189 : real(dp), contiguous, pointer :: cwavef(:,:,:) => null()
190 : real(dp), contiguous, pointer :: cwavefb(:,:,:) => null()
191 : real(dp), contiguous, pointer :: cwavef_x(:,:) => null()
192 : real(dp), contiguous, pointer :: cwavef_y(:,:) => null()
193 : real(dp), contiguous, pointer :: cwavefb_2(:,:) => null() ! only use when paral_kgb=0
194 : real(dp), contiguous, pointer :: cwavefb_x(:,:) => null() ! only use when paral_kgb=0
195 : real(dp), contiguous, pointer :: cwavefb_y(:,:) => null() ! only use when paral_kgb=0
196 : real(dp), contiguous, pointer :: rhoaug(:,:,:) => null()
197 : real(dp), contiguous, pointer :: rhoaug_down(:,:,:) => null()
198 : real(dp), contiguous, pointer :: rhoaug_up(:,:,:) => null()
199 : real(dp), contiguous, pointer :: rhoaug_mx(:,:,:) => null()
200 : real(dp), contiguous, pointer :: rhoaug_my(:,:,:) => null()
201 : real(dp), contiguous, pointer :: wfraug(:,:,:,:) => null()
202 : real(dp), contiguous, pointer :: cg_k(:,:) => null()
203 : ! *************************************************************************
204 :
205 : DBG_ENTER("COLL")
206 :
207 17150 : call timab(790+tim_mkrho,1,tsec)
208 17150 : call timab(799,1,tsec)
209 :
210 17150 : if(mpi_enreg%paralbd==0) tim_fourwf=3
211 17150 : if(mpi_enreg%paralbd==1) tim_fourwf=6
212 :
213 17150 : if(.not.(present(option))) then
214 16308 : ioption=0
215 : else
216 842 : ioption=option
217 : end if
218 :
219 17150 : if(.not.(present(printout))) then
220 : l_printout = .true.
221 : else
222 1205 : l_printout = printout
223 : end if
224 :
225 : ! Not sure what to do for Wannier90 DMFT
226 17150 : if(ioption/=0.and.(paw_dmft%use_sc_dmft==1.or.paw_dmft%use_sc_dmft==10)) then
227 0 : ABI_ERROR('option argument value of this routines should be 0 if usedmft=1 or 10.')
228 : end if
229 17150 : if(paw_dmft%use_sc_dmft/=0.and.paw_dmft%use_sc_dmft/=10) then
230 138 : nbandc1=(paw_dmft%mbandc-1)*paw_dmft%use_sc_dmft+1
231 17012 : else if(paw_dmft%use_sc_dmft==10) then
232 0 : nbandc1=paw_dmft%mbandc
233 : else
234 : nbandc1=1
235 : end if
236 17150 : use_nondiag_occup_dmft=0
237 :
238 : !if(dtset%nspinor==2.and.paw_dmft%use_sc_dmft==1) then
239 : !write(message, '(a,a,a,a)' )ch10,&
240 : !& ' mkrho : ERROR -',ch10,&
241 : !& ' nspinor argument value of this routines should be 1 if usedmft=1. '
242 : !call wrtout(std_out,message,'COLL')
243 : !end if
244 :
245 17150 : my_nspinor=max(1,dtset%nspinor/mpi_enreg%nproc_spinor)
246 17150 : if (mpi_enreg%paral_spinor==0) then
247 17146 : ispinor_index=1;jspinor_index=1
248 17146 : nspinor1TreatedByThisProc=.true.
249 17146 : nspinor2TreatedByThisProc=(dtset%nspinor==2)
250 : else
251 4 : ispinor_index=mpi_enreg%me_spinor+1;jspinor_index=3-ispinor_index
252 4 : nspinor1TreatedByThisProc=(mpi_enreg%me_spinor==0)
253 4 : nspinor2TreatedByThisProc=(mpi_enreg%me_spinor==1)
254 : end if
255 :
256 : !Set local variable which depend on option argument
257 :
258 : !nalpha*nbeta is the number of element of the kinetic energy density tensor
259 : !to be computed in the irreducible Brillouin Zone (BZ) to get the result in the full BZ.
260 : !In case of electron density calculation, nalpha=nbeta=1
261 : select case (ioption)
262 : case (0)
263 : nalpha = 1
264 824 : nbeta = 1
265 : case (1)
266 824 : nalpha = 3
267 824 : nbeta = 1
268 4120 : ABI_MALLOC(taur_alphabeta,(dtset%nfft,dtset%nspden,3,1))
269 : case (2)
270 0 : nalpha = 3
271 0 : nbeta = 3
272 0 : ABI_MALLOC(taur_alphabeta,(dtset%nfft,dtset%nspden,3,3))
273 : case default
274 17150 : ABI_BUG(sjoin('ioption argument value should be 0,1 or 2 while got:', itoa(ioption)))
275 : end select
276 :
277 : !Init me
278 17150 : me=mpi_enreg%me_kpt
279 : !zero the charge density array in real space
280 : !$OMP PARALLEL DO COLLAPSE(2)
281 42053 : do ispden=1,dtset%nspden
282 232562881 : do ifft=1,dtset%nfft
283 232545731 : rhor(ifft,ispden)=zero
284 : end do
285 : end do
286 :
287 : !WVL - Branching with a separate mkrho procedure in wavelet.
288 17150 : if (dtset%usewvl == 1) then
289 0 : select case(ioption)
290 : case (0)
291 0 : call wvl_mkrho(dtset, irrzon, mpi_enreg, phnons, rhor, wvl_wfs, wvl_den)
292 0 : return
293 : case (1)
294 : !call wvl_mkrho(dtset, mpi_enreg, occ, rhor, wvl_wfs, wvl_den)
295 0 : ABI_ERROR("kinetic energy density (taur) is not yet implemented in wavelet formalism.")
296 : case (2)
297 : !call wvl_mkrho(dtset, mpi_enreg, occ, rhor, wvl_wfs, wvl_den)
298 0 : ABI_BUG('kinetic energy density tensor (taur_(alpha,beta)) is not yet implemented in wavelet formalism.')
299 : end select
300 : end if
301 : !WVL - Following is done in plane waves.
302 :
303 17150 : gpu_option=dtset%gpu_option
304 17150 : gpu_cwavef=(gpu_option==ABI_GPU_OPENMP .and. paw_dmft%use_sc_dmft/=1)
305 17150 : transfer_cg = .false.
306 : #ifdef HAVE_OPENMP_OFFLOAD
307 : if(gpu_option==ABI_GPU_OPENMP) then
308 : transfer_cg = .not. xomp_target_is_present(c_loc(cg))
309 : if(mpi_enreg%paral_kgb==0) then
310 : !FIXME If cg is on card, we need to rapatriate it on CPU when paral_kgb==0
311 : !$OMP TARGET UPDATE FROM(cg) IF(.not. transfer_cg)
312 : end if
313 : end if
314 : #endif
315 :
316 17150 : l_nfft_blocks=1; if(present(nfft_blocks) .and. gpu_option/=ABI_GPU_DISABLED) l_nfft_blocks=nfft_blocks
317 :
318 : !start loop over alpha and beta
319 :
320 35948 : do alpha=1,nalpha
321 54746 : do beta=1,nbeta
322 :
323 : ! start loop over spin and k points
324 18798 : bdtot_index=0
325 18798 : icg=0
326 :
327 : ! n4,n5,n6 are FFT dimensions, modified to avoid cache trashing
328 18798 : n1 = dtset%ngfft(1) ; n2 = dtset%ngfft(2) ; n3 = dtset%ngfft(3)
329 18798 : n4 = dtset%ngfft(4) ; n5 = dtset%ngfft(5) ; n6 = dtset%ngfft(6)
330 :
331 18798 : ndat = 1
332 18798 : if (mpi_enreg%paral_kgb==1) then
333 4488 : ndat = mpi_enreg%bandpp
334 14310 : else if (gpu_option/=ABI_GPU_DISABLED) then
335 0 : ndat = dtset%mband
336 : end if
337 :
338 18798 : if (gpu_option == ABI_GPU_KOKKOS) then
339 : #if defined HAVE_GPU && defined HAVE_YAKL
340 : ABI_MALLOC_MANAGED(cwavef,(/2,dtset%mpw*ndat,my_nspinor/))
341 : #endif
342 : else
343 18798 : if (gpu_option/=ABI_GPU_DISABLED) then
344 0 : ABI_MALLOC(cwavef,(2,dtset%mpw*ndat,my_nspinor))
345 : else
346 75192 : ABI_MALLOC(cwavef,(2,dtset%mpw,my_nspinor))
347 : end if
348 : end if
349 :
350 18798 : if(gpu_option == ABI_GPU_KOKKOS) then
351 : #if defined HAVE_GPU && defined HAVE_YAKL
352 : ABI_MALLOC_MANAGED(rhoaug, (/n4,n5,n6/))
353 : ABI_MALLOC_MANAGED(wfraug, (/2,n4,n5,n6*ndat/))
354 : ABI_MALLOC_MANAGED(cwavefb, (/2,dtset%mpw*paw_dmft%use_sc_dmft,my_nspinor/))
355 : if(dtset%nspden==4) then
356 : ABI_MALLOC_MANAGED(rhoaug_up, (/n4,n5,n6/))
357 : ABI_MALLOC_MANAGED(rhoaug_down,(/n4,n5,n6/))
358 : ABI_MALLOC_MANAGED(rhoaug_mx, (/n4,n5,n6/))
359 : ABI_MALLOC_MANAGED(rhoaug_my, (/n4,n5,n6/))
360 : rhoaug_up(:,:,:)=zero
361 : rhoaug_down(:,:,:)=zero
362 : rhoaug_mx(:,:,:)=zero
363 : rhoaug_my(:,:,:)=zero
364 : end if
365 : #endif
366 : else
367 93990 : ABI_MALLOC(rhoaug, (n4,n5,n6))
368 : #ifdef HAVE_OPENMP_OFFLOAD
369 : !$OMP TARGET ENTER DATA MAP(alloc:rhoaug) IF(gpu_option==ABI_GPU_OPENMP)
370 : #endif
371 93990 : ABI_MALLOC(wfraug, (2,n4,n5,n6*(ndat/l_nfft_blocks+ndat-(ndat/l_nfft_blocks)*l_nfft_blocks)))
372 75192 : ABI_MALLOC(cwavefb, (2,dtset%mpw*paw_dmft%use_sc_dmft,my_nspinor))
373 18798 : if(dtset%nspden==4) then
374 7828 : ABI_MALLOC(rhoaug_up, (n4,n5,n6))
375 7828 : ABI_MALLOC(rhoaug_down,(n4,n5,n6))
376 7828 : ABI_MALLOC(rhoaug_mx, (n4,n5,n6))
377 7828 : ABI_MALLOC(rhoaug_my, (n4,n5,n6))
378 14522667 : rhoaug_up(:,:,:)=zero
379 14522667 : rhoaug_down(:,:,:)=zero
380 14522667 : rhoaug_mx(:,:,:)=zero
381 14522667 : rhoaug_my(:,:,:)=zero
382 : end if
383 : end if
384 :
385 39779 : do isppol=1,dtset%nsppol
386 20981 : ikg=0
387 :
388 20981 : if(gpu_option==ABI_GPU_OPENMP) then
389 0 : call gpu_set_to_zero(rhoaug,int(n4,c_size_t)*n5*n6)
390 : else
391 305532680 : rhoaug(:,:,:)=zero
392 : end if
393 240846 : do ikpt=1,dtset%nkpt
394 :
395 219865 : nband_k = dtset%nband(ikpt+(isppol-1)*dtset%nkpt)
396 219865 : mband_mem = nband_k
397 219865 : if (dtset%paral_kgb==0) mband_mem = nband_k/mpi_enreg%nproc_band
398 219865 : npw_k=npwarr(ikpt)
399 219865 : istwf_k = dtset%istwfk(ikpt)
400 :
401 219865 : if(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,isppol,me)) then
402 46171 : bdtot_index=bdtot_index+nband_k
403 46171 : cycle
404 : end if
405 :
406 694776 : ABI_MALLOC(gbound,(2*dtset%mgfft+8,2))
407 173694 : if(gpu_option == ABI_GPU_KOKKOS) then
408 : #if defined HAVE_GPU && defined HAVE_YAKL
409 : ABI_MALLOC_MANAGED(kg_k, (/3,npw_k/))
410 : #endif
411 : else
412 521082 : ABI_MALLOC(kg_k,(3,npw_k))
413 : end if
414 :
415 144775570 : kg_k(:,1:npw_k)=kg(:,1+ikg:npw_k+ikg)
416 173694 : call sphereboundary(gbound,istwf_k,kg_k,dtset%mgfft,npw_k)
417 :
418 : ! Loop over bands to fft and square for rho(r)
419 : ! Should be changed to treat bands by batch always
420 :
421 173694 : if (mpi_enreg%paral_kgb /= 1) then ! Not yet parallelized on spinors
422 :
423 : #ifdef HAVE_OPENMP_OFFLOAD
424 : ! With OpenMP GPU, uploading kg_k when paral_kgb==0
425 : !$OMP TARGET ENTER DATA MAP(to:kg_k) IF(gpu_option==ABI_GPU_OPENMP)
426 : #endif
427 :
428 157686 : if (gpu_option /= ABI_GPU_DISABLED .and. paw_dmft%use_sc_dmft/=1) then
429 : !On GPU, treat all bands at once
430 0 : ABI_MALLOC(weight_t,(nband_k))
431 0 : nband_occ = 0
432 0 : do iband=1,nband_k
433 0 : ipwsp = (iband-1)*npw_k*my_nspinor + icg
434 0 : locc_test = abs(occ(iband+bdtot_index))>tol8
435 0 : if (locc_test) then
436 0 : nband_occ = nband_occ +1
437 0 : ipwbd = (nband_occ-1) * npw_k
438 0 : weight_t(nband_occ) = occ(iband+bdtot_index) * dtset%wtk(ikpt)/ucvol
439 0 : cwavef(:,ipwbd+1:ipwbd+npw_k,1) = cg(:,ipwsp+1:ipwsp+npw_k)
440 0 : if (my_nspinor==2) cwavef(:,ipwbd+1:ipwbd+npw_k,2) = cg(:,ipwsp+npw_k+1:ipwsp+npw_k+npw_k)
441 0 : if (ioption==1) then ! Multiplication by 2pi i (k+G)_alpha
442 0 : gp2pi1 = gprimd(alpha,1)*two_pi ; gp2pi2 = gprimd(alpha,2)*two_pi ; gp2pi3 = gprimd(alpha,3)*two_pi
443 0 : kpt_cart = gp2pi1*dtset%kptns(1,ikpt) + gp2pi2*dtset%kptns(2,ikpt) + gp2pi3*dtset%kptns(3,ikpt)
444 0 : do ispinor=1,my_nspinor
445 0 : do ipw=1,npw_k
446 0 : kg_k_cart = gp2pi1*kg_k(1,ipw) + gp2pi2*kg_k(2,ipw) + gp2pi3*kg_k(3,ipw) + kpt_cart
447 0 : cwftmp = -cwavef(2,ipwbd+ipw,ispinor)*kg_k_cart
448 0 : cwavef(2,ipwbd+ipw,ispinor) = cwavef(1,ipwbd+ipw,ispinor)*kg_k_cart
449 0 : cwavef(1,ipwbd+ipw,ispinor) = cwftmp
450 : end do
451 : end do
452 0 : else if (ioption==2) then
453 0 : ABI_ERROR('kinetic energy density tensor (taur_(alpha,beta)) is not yet implemented.')
454 : end if ! end if ioption==1
455 : end if ! end if locc_test
456 : end do ! end iband=1,nband_k
457 0 : if (nband_occ>0) then
458 0 : nfft_blocks_occ=1; nband_fftblock_occ=nband_occ
459 0 : if(l_nfft_blocks>1) then
460 0 : nfft_blocks_occ = real(nband_occ)/mpi_enreg%bandpp * l_nfft_blocks
461 0 : nfft_blocks_occ = min(nband_occ,nfft_blocks_occ)
462 0 : nband_fftblock_occ = nband_occ/nfft_blocks_occ
463 : end if
464 :
465 : call fourwf_optmem(1,rhoaug,cwavef(:,1:nband_occ*npw_k,1),dummy,wfraug(:,:,:,1:n6*nband_fftblock_occ),&
466 : & gbound,gbound,istwf_k,kg_k,kg_k,dtset%mgfft,mpi_enreg,nband_occ,nfft_blocks_occ,dtset%ngfft,&
467 : & npw_k,1,n4,n5,n6,1,tim_fourwf,weight,weight_i,&
468 : & weight_array_r=weight_t(1:nband_occ),weight_array_i=weight_t(1:nband_occ),&
469 0 : & gpu_option=gpu_option)
470 0 : if(dtset%nspinor==2)then
471 0 : if(dtset%nspden==1) then
472 : ! We need only the total density : accumulation continues on top of rhoaug
473 : call fourwf_optmem(1,rhoaug,cwavef(:,1:nband_occ*npw_k,2),dummy,wfraug(:,:,:,1:n6*nband_fftblock_occ),&
474 : & gbound,gbound,istwf_k,kg_k,kg_k,dtset%mgfft,mpi_enreg,nband_occ,nfft_blocks_occ,dtset%ngfft,&
475 : & npw_k,1,n4,n5,n6,1,tim_fourwf,weight,weight_i,&
476 : & weight_array_r=weight_t(1:nband_occ),weight_array_i=weight_t(1:nband_occ),&
477 0 : & gpu_option=gpu_option)
478 0 : else if(dtset%nspden==4) then
479 : ! Build the four components of rho. We use only norm quantities and, so fourwf.
480 : ! $\sum_{n} f_n \Psi^{* \alpha}_n \Psi^{\alpha}_n =\rho^{\alpha \alpha}$
481 : ! $\sum_{n} f_n (\Psi^{1}+\Psi^{2})^*_n (\Psi^{1}+\Psi^{2})_n=rho+m_x$
482 : ! $\sum_{n} f_n (\Psi^{1}-i \Psi^{2})^*_n (\Psi^{1}-i \Psi^{2})_n=rho+m_y$
483 0 : if(gpu_option == ABI_GPU_KOKKOS) then
484 : #if defined HAVE_GPU && defined HAVE_YAKL
485 : ABI_MALLOC_MANAGED(cwavef_x, (/2,npw_k*nband_occ/))
486 : ABI_MALLOC_MANAGED(cwavef_y, (/2,npw_k*nband_occ/))
487 : ABI_MALLOC_MANAGED(cwavefb_2,(/2,npw_k*nband_occ*paw_dmft%use_sc_dmft/))
488 : ABI_MALLOC_MANAGED(cwavefb_x,(/2,npw_k*nband_occ*paw_dmft%use_sc_dmft/))
489 : ABI_MALLOC_MANAGED(cwavefb_y,(/2,npw_k*nband_occ*paw_dmft%use_sc_dmft/))
490 : #endif
491 : else
492 0 : ABI_MALLOC(cwavef_x,(2,npw_k*nband_occ))
493 0 : ABI_MALLOC(cwavef_y,(2,npw_k*nband_occ))
494 0 : ABI_MALLOC(cwavefb_2,(2,npw_k*nband_occ*paw_dmft%use_sc_dmft))
495 0 : ABI_MALLOC(cwavefb_x,(2,npw_k*nband_occ*paw_dmft%use_sc_dmft))
496 0 : ABI_MALLOC(cwavefb_y,(2,npw_k*nband_occ*paw_dmft%use_sc_dmft))
497 : end if
498 : ! $(\Psi^{1}+\Psi^{2})$
499 0 : cwavef_x(:,:)=cwavef(:,1:npw_k*nband_occ,1)+cwavef(:,1:npw_k*nband_occ,2)
500 : ! $(\Psi^{1}-i \Psi^{2})$
501 0 : cwavef_y(1,:)=cwavef(1,1:npw_k*nband_occ,1)+cwavef(2,1:npw_k*nband_occ,2)
502 0 : cwavef_y(2,:)=cwavef(2,1:npw_k*nband_occ,1)-cwavef(1,1:npw_k*nband_occ,2)
503 0 : if(use_nondiag_occup_dmft==1) then
504 0 : cwavefb_2(:,:)=cwavefb(:,1:npw_k*nband_occ,2)
505 0 : cwavefb_x(:,:)=cwavefb(:,1:npw_k*nband_occ,1)+cwavefb(:,1:npw_k*nband_occ,2)
506 0 : cwavefb_y(1,:)=cwavefb(1,1:npw_k*nband_occ,1)+cwavefb(2,1:npw_k*nband_occ,2)
507 0 : cwavefb_y(2,:)=cwavefb(2,1:npw_k*nband_occ,1)-cwavefb(1,1:npw_k*nband_occ,2)
508 : end if
509 : #ifdef HAVE_OPENMP_OFFLOAD
510 : !$OMP TARGET UPDATE FROM(rhoaug) IF(gpu_option==ABI_GPU_OPENMP)
511 : #endif
512 0 : rhoaug_up(:,:,:)=rhoaug(:,:,:) !Already computed
513 : call fourwf_optmem(1,rhoaug_down,cwavef(:,1:nband_occ*npw_k,2),dummy,wfraug(:,:,:,1:n6*nband_fftblock_occ),gbound,gbound,&
514 : & istwf_k,kg_k,kg_k,dtset%mgfft,mpi_enreg,nband_occ,nfft_blocks_occ,dtset%ngfft,&
515 : & npw_k,1,n4,n5,n6,1,tim_fourwf,weight,weight_i,&
516 : & weight_array_r=weight_t(1:nband_occ),weight_array_i=weight_t(1:nband_occ),&
517 : & use_ndo=use_nondiag_occup_dmft,fofginb=cwavefb_2,&
518 0 : & gpu_option=gpu_option)
519 :
520 : call fourwf_optmem(1,rhoaug_mx,cwavef_x,dummy,wfraug(:,:,:,1:n6*nband_fftblock_occ),gbound,gbound,&
521 : & istwf_k,kg_k,kg_k,dtset%mgfft,mpi_enreg,nband_occ,nfft_blocks_occ,dtset%ngfft,&
522 : & npw_k,1,n4,n5,n6,1,tim_fourwf,weight,weight_i,&
523 : & weight_array_r=weight_t(1:nband_occ),weight_array_i=weight_t(1:nband_occ),&
524 : & use_ndo=use_nondiag_occup_dmft,fofginb=cwavefb_x,&
525 0 : & gpu_option=gpu_option)
526 :
527 : call fourwf_optmem(1,rhoaug_my,cwavef_y,dummy,wfraug(:,:,:,1:n6*nband_fftblock_occ),gbound,gbound,&
528 : & istwf_k,kg_k,kg_k,dtset%mgfft,mpi_enreg,nband_occ,nfft_blocks_occ,dtset%ngfft,&
529 : & npw_k,1,n4,n5,n6,1,tim_fourwf,weight,weight_i,&
530 : & weight_array_r=weight_t(1:nband_occ),weight_array_i=weight_t(1:nband_occ),&
531 : & use_ndo=use_nondiag_occup_dmft,fofginb=cwavefb_y,&
532 0 : & gpu_option=gpu_option)
533 0 : if(gpu_option == ABI_GPU_KOKKOS) then
534 : #if defined HAVE_GPU && defined HAVE_YAKL
535 : ABI_FREE_MANAGED(cwavef_x)
536 : ABI_FREE_MANAGED(cwavef_y)
537 : ABI_FREE_MANAGED(cwavefb_2)
538 : ABI_FREE_MANAGED(cwavefb_x)
539 : ABI_FREE_MANAGED(cwavefb_y)
540 : #endif
541 : else
542 0 : ABI_FREE(cwavef_x)
543 0 : ABI_FREE(cwavef_y)
544 0 : ABI_FREE(cwavefb_2)
545 0 : ABI_FREE(cwavefb_x)
546 0 : ABI_FREE(cwavefb_y)
547 : end if
548 : end if
549 : end if
550 : end if
551 0 : ABI_FREE(weight_t)
552 :
553 : else ! CPU version
554 :
555 157686 : iband_me = 0
556 1649267 : do iband=1,nband_k
557 1491581 : if(mpi_enreg%paralbd==1)then
558 259143 : if(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,iband,iband,isppol,me)) cycle
559 : end if
560 1491437 : iband_me = iband_me + 1
561 3332284 : do ibandc1=1,nbandc1 ! in case of DMFT
562 : ! Check if DMFT and only treat occupied states (check on occ.)
563 1683161 : if(paw_dmft%use_sc_dmft == 1 .or. paw_dmft%use_sc_dmft == 10) then
564 206950 : iband1 = paw_dmft%include_bands(ibandc1)
565 206950 : if(paw_dmft%band_in(iband)) then
566 164526 : if(.not. paw_dmft%band_in(iband1)) stop
567 164526 : use_nondiag_occup_dmft = 1
568 : locc_test = abs(paw_dmft%occnd(1,iband,iband1,ikpt,isppol)) +&
569 164526 : & abs(paw_dmft%occnd(2,iband,iband1,ikpt,isppol))>tol8
570 : else
571 42424 : use_nondiag_occup_dmft = 0
572 42424 : locc_test = abs(occ(iband+bdtot_index))>tol8
573 42424 : if(ibandc1 /=1 .and. .not. paw_dmft%band_in(iband)) cycle
574 : end if
575 : else
576 1476211 : use_nondiag_occup_dmft = 0
577 1476211 : locc_test = abs(occ(iband+bdtot_index))>tol8
578 : end if
579 3139398 : if (locc_test) then
580 : ! Obtain Fourier transform in fft box and accumulate the density or
581 : ! the kinetic energy density
582 : ! Not yet parallelized on nspinor if paral_kgb/=1
583 964171 : ipwsp=(iband_me-1)*npw_k*my_nspinor +icg
584 1576517569 : cwavef(:,1:npw_k,1) = cg(:,1+ipwsp :ipwsp+npw_k)
585 223374745 : if (my_nspinor==2) cwavef(:,1:npw_k,2)=cg(:,1+ipwsp+npw_k:ipwsp+2*npw_k)
586 964171 : if(ioption==1)then
587 : ! Multiplication by 2pi i (k+G)_alpha
588 26577 : gp2pi1=gprimd(alpha,1)*two_pi ; gp2pi2=gprimd(alpha,2)*two_pi ; gp2pi3=gprimd(alpha,3)*two_pi
589 26577 : kpt_cart=gp2pi1*dtset%kptns(1,ikpt)+gp2pi2*dtset%kptns(2,ikpt)+gp2pi3*dtset%kptns(3,ikpt)
590 63978 : do ispinor=1,my_nspinor
591 19228755 : do ipw=1,npw_k
592 19164777 : kg_k_cart=gp2pi1*kg_k(1,ipw)+gp2pi2*kg_k(2,ipw)+gp2pi3*kg_k(3,ipw)+kpt_cart
593 19164777 : cwftmp=-cwavef(2,ipw,ispinor)*kg_k_cart
594 19164777 : cwavef(2,ipw,ispinor)=cwavef(1,ipw,ispinor)*kg_k_cart
595 19202178 : cwavef(1,ipw,ispinor)=cwftmp
596 : end do
597 : end do
598 937594 : else if(ioption==2)then
599 0 : ABI_ERROR('kinetic energy density tensor (taur_(alpha,beta)) is not yet implemented.')
600 : end if
601 : ! Non diag occupation in DMFT.
602 : ! TODO : this will break in full distrib of band memory
603 964171 : if(use_nondiag_occup_dmft==1) then
604 28108 : ipwsp=(iband1-1)*npw_k*my_nspinor +icg
605 59241976 : cwavefb(:,1:npw_k,1)=cg(:,1+ipwsp:ipwsp+npw_k)
606 2009404 : if (my_nspinor==2) cwavefb(:,1:npw_k,2)=cg(:,ipwsp+npw_k+1:ipwsp+2*npw_k)
607 28108 : weight =paw_dmft%occnd(1,iband,iband1,ikpt,isppol)*dtset%wtk(ikpt)/ucvol
608 28108 : weight_i=paw_dmft%occnd(2,iband,iband1,ikpt,isppol)*dtset%wtk(ikpt)/ucvol
609 : else
610 936063 : weight=occ(iband+bdtot_index)*dtset%wtk(ikpt)/ucvol
611 936063 : weight_i=weight
612 : end if
613 :
614 : ! The same section of code is also found in vtowfk.F90 : should be rationalized !
615 :
616 : call fourwf(1,rhoaug,cwavef(:,:,1),dummy,wfraug,gbound,gbound,&
617 : & istwf_k,kg_k,kg_k,dtset%mgfft,mpi_enreg,1,dtset%ngfft,&
618 : & npw_k,1,n4,n5,n6,1,tim_fourwf,weight,weight_i,&
619 : & use_ndo=use_nondiag_occup_dmft,fofginb=cwavefb(:,:,1),&
620 964171 : & gpu_option=gpu_option)
621 964171 : if(dtset%nspinor==2)then
622 172481 : if(dtset%nspden==1) then
623 : ! We need only the total density : accumulation continues on top of rhoaug
624 : call fourwf(1,rhoaug,cwavef(:,:,2),dummy,wfraug,gbound,gbound,&
625 : & istwf_k,kg_k,kg_k,dtset%mgfft,mpi_enreg,1,dtset%ngfft,&
626 : & npw_k,1,n4,n5,n6,1,tim_fourwf,weight,weight_i,&
627 : & use_ndo=use_nondiag_occup_dmft,fofginb=cwavefb(:,:,2),&
628 41164 : & gpu_option=gpu_option)
629 131317 : else if(dtset%nspden==4) then
630 : ! Build the four components of rho. We use only norm quantities and, so fourwf.
631 : ! $\sum_{n} f_n \Psi^{* \alpha}_n \Psi^{\alpha}_n =\rho^{\alpha \alpha}$
632 : ! $\sum_{n} f_n (\Psi^{1}+\Psi^{2})^*_n (\Psi^{1}+\Psi^{2})_n=rho+m_x$
633 : ! $\sum_{n} f_n (\Psi^{1}-i \Psi^{2})^*_n (\Psi^{1}-i \Psi^{2})_n=rho+m_y$
634 131317 : if(gpu_option == ABI_GPU_KOKKOS) then
635 : #if defined HAVE_GPU && defined HAVE_YAKL
636 : ABI_MALLOC_MANAGED(cwavef_x, (/2,npw_k/))
637 : ABI_MALLOC_MANAGED(cwavef_y, (/2,npw_k/))
638 : ABI_MALLOC_MANAGED(cwavefb_x,(/2,npw_k*paw_dmft%use_sc_dmft/))
639 : ABI_MALLOC_MANAGED(cwavefb_y,(/2,npw_k*paw_dmft%use_sc_dmft/))
640 : #endif
641 : else
642 393951 : ABI_MALLOC(cwavef_x,(2,npw_k))
643 262634 : ABI_MALLOC(cwavef_y,(2,npw_k))
644 393951 : ABI_MALLOC(cwavefb_x,(2,npw_k*paw_dmft%use_sc_dmft))
645 262634 : ABI_MALLOC(cwavefb_y,(2,npw_k*paw_dmft%use_sc_dmft))
646 : end if
647 : ! $(\Psi^{1}+\Psi^{2})$
648 155960035 : cwavef_x(:,:)=cwavef(:,1:npw_k,1)+cwavef(:,1:npw_k,2)
649 : ! $(\Psi^{1}-i \Psi^{2})$
650 52074223 : cwavef_y(1,:)=cwavef(1,1:npw_k,1)+cwavef(2,1:npw_k,2)
651 52074223 : cwavef_y(2,:)=cwavef(2,1:npw_k,1)-cwavef(1,1:npw_k,2)
652 131317 : if(use_nondiag_occup_dmft==1) then
653 1982930 : cwavefb_x(:,:)=cwavefb(:,1:npw_k,1)+cwavefb(:,1:npw_k,2)
654 662066 : cwavefb_y(1,:)=cwavefb(1,1:npw_k,1)+cwavefb(2,1:npw_k,2)
655 662066 : cwavefb_y(2,:)=cwavefb(2,1:npw_k,1)-cwavefb(1,1:npw_k,2)
656 : end if
657 : #ifdef HAVE_OPENMP_OFFLOAD
658 : !$OMP TARGET UPDATE FROM(rhoaug) IF(gpu_option==ABI_GPU_OPENMP)
659 : #endif
660 1460911819 : rhoaug_up(:,:,:)=rhoaug(:,:,:) !Already computed
661 : call fourwf(1,rhoaug_down,cwavef(:,:,2),dummy,wfraug,gbound,gbound,&
662 : & istwf_k,kg_k,kg_k,dtset%mgfft,mpi_enreg,1,dtset%ngfft,&
663 : & npw_k,1,n4,n5,n6,1,tim_fourwf,weight,weight_i,&
664 : & use_ndo=use_nondiag_occup_dmft,fofginb=cwavefb(:,:,2),&
665 131317 : & gpu_option=gpu_option)
666 :
667 : call fourwf(1,rhoaug_mx,cwavef_x,dummy,wfraug,gbound,gbound,&
668 : & istwf_k,kg_k,kg_k,dtset%mgfft,mpi_enreg,1,dtset%ngfft,&
669 : & npw_k,1,n4,n5,n6,1,tim_fourwf,weight,weight_i,&
670 : & use_ndo=use_nondiag_occup_dmft,fofginb=cwavefb_x,&
671 131317 : & gpu_option=gpu_option)
672 :
673 : call fourwf(1,rhoaug_my,cwavef_y,dummy,wfraug,gbound,gbound,&
674 : & istwf_k,kg_k,kg_k,dtset%mgfft,mpi_enreg,1,dtset%ngfft,&
675 : & npw_k,1,n4,n5,n6,1,tim_fourwf,weight,weight_i,&
676 : & use_ndo=use_nondiag_occup_dmft,fofginb=cwavefb_y,&
677 131317 : & gpu_option=gpu_option)
678 131317 : if(gpu_option == ABI_GPU_KOKKOS) then
679 : #if defined HAVE_GPU && defined HAVE_YAKL
680 : ABI_FREE_MANAGED(cwavef_x)
681 : ABI_FREE_MANAGED(cwavef_y)
682 : ABI_FREE_MANAGED(cwavefb_x)
683 : ABI_FREE_MANAGED(cwavefb_y)
684 : #endif
685 : else
686 131317 : ABI_FREE(cwavef_x)
687 131317 : ABI_FREE(cwavef_y)
688 131317 : ABI_FREE(cwavefb_x)
689 131317 : ABI_FREE(cwavefb_y)
690 : end if
691 : end if ! dtset%nspden/=4
692 : end if
693 : else
694 : ! Accumulate the number of one-way 3D ffts skipped
695 683646 : nskip=nskip+1
696 : end if ! abs(occ(iband+bdtot_index))>tol8
697 : end do ! iband1=1,(nband_k-1)*paw_dmft%use_sc_dmft+1
698 : end do ! iband=1,nband_k
699 :
700 : end if ! gpu_option
701 :
702 :
703 : #ifdef HAVE_OPENMP_OFFLOAD
704 : !$OMP TARGET EXIT DATA MAP(delete:kg_k) IF(gpu_option==ABI_GPU_OPENMP)
705 : #endif
706 : else !paral_kgb==1
707 :
708 16008 : if(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,isppol,me)) cycle
709 :
710 16008 : call bandfft_kpt_set_ikpt(ikpt,mpi_enreg)
711 16008 : nbdblock=nband_k/(mpi_enreg%nproc_band * mpi_enreg%bandpp)
712 16008 : blocksize=nband_k/nbdblock
713 16008 : cg_k => cg(:,1+icg:npw_k*my_nspinor*blocksize*nbdblock+icg)
714 : #ifdef HAVE_OPENMP_OFFLOAD
715 : !$OMP TARGET ENTER DATA MAP(to:cg_k) if(gpu_cwavef .and. transfer_cg)
716 : #endif
717 :
718 16008 : if(gpu_option == ABI_GPU_KOKKOS) then
719 : #if defined HAVE_GPU && defined HAVE_YAKL
720 : if(associated(cwavef)) then
721 : ABI_FREE_MANAGED(cwavef)
722 : end if
723 : ABI_MALLOC_MANAGED(cwavef,(/2,npw_k*blocksize,dtset%nspinor/))
724 : #endif
725 : else
726 16008 : if(associated(cwavef)) then
727 4488 : ABI_FREE(cwavef)
728 : end if
729 64032 : ABI_MALLOC(cwavef,(2,npw_k*blocksize,dtset%nspinor))
730 : #ifdef HAVE_OPENMP_OFFLOAD
731 : !$OMP TARGET ENTER DATA MAP(alloc:cwavef) IF(gpu_cwavef)
732 : #endif
733 : end if
734 16008 : if(ioption==1) then
735 3618 : ABI_MALLOC(kg_k_cart_block,(npw_k))
736 : end if
737 48024 : ABI_MALLOC(occ_k,(nband_k))
738 241128 : occ_k(:)=occ(bdtot_index+1:bdtot_index+nband_k)
739 :
740 : ! ---------- DMFT
741 16008 : if(allocated(cwavef_rot)) then
742 176 : ABI_FREE(cwavef_rot)
743 176 : ABI_FREE(occ_diag)
744 : ! ABI_FREE(occ_nd)
745 : end if
746 16008 : if(paw_dmft%use_sc_dmft==1.or.paw_dmft%use_sc_dmft==10) then
747 : ! Allocation of DMFT temporaries arrays
748 1100 : ABI_MALLOC(cwavef_rot,(2,npw_k,blocksize,dtset%nspinor))
749 660 : ABI_MALLOC(occ_diag,(blocksize))
750 : ! ABI_MALLOC(occ_nd,(2, blocksize, blocksize, dtset%nspinor))
751 : end if
752 : ! ---------- END DMFT
753 :
754 95240 : do iblock=1,nbdblock
755 79232 : if (dtset%nspinor==1) then
756 73576 : if(gpu_cwavef) then
757 : #ifdef HAVE_OPENMP_OFFLOAD
758 : call gpu_copy(cwavef(:,1:npw_k*blocksize,1), &
759 : & cg_k(:,1+(iblock-1)*npw_k*blocksize:iblock*npw_k*blocksize), &
760 : & int(2,c_size_t)*npw_k*blocksize)
761 : #endif
762 : else
763 174637072 : cwavef(:,1:npw_k*blocksize,1)=cg_k(:,1+(iblock-1)*npw_k*blocksize:iblock*npw_k*blocksize)
764 : end if
765 : else
766 5656 : if (mpi_enreg%paral_spinor==0) then
767 4856 : ishf=(iblock-1)*npw_k*my_nspinor*blocksize
768 4856 : if(gpu_cwavef) then
769 : #ifdef HAVE_OPENMP_OFFLOAD
770 : !$OMP TARGET TEAMS DISTRIBUTE MAP(to:cg_k,cwavef) PRIVATE(ib)
771 : do ib=1,blocksize
772 : !$OMP PARALLEL DO PRIVATE(ipw)
773 : do ipw=1,npw_k
774 : cwavef(1,ipw+(ib-1)*npw_k,1)=cg_k(1,ipw+(2*ib-2)*npw_k+ishf)
775 : cwavef(2,ipw+(ib-1)*npw_k,1)=cg_k(2,ipw+(2*ib-2)*npw_k+ishf)
776 :
777 : cwavef(1,ipw+(ib-1)*npw_k,2)=cg_k(1,ipw+(2*ib-1)*npw_k+ishf)
778 : cwavef(2,ipw+(ib-1)*npw_k,2)=cg_k(2,ipw+(2*ib-1)*npw_k+ishf)
779 : end do
780 : end do
781 : #endif
782 : else
783 26296 : do ib=1,blocksize
784 27266240 : cwavef(:,(ib-1)*npw_k+1:ib*npw_k,1)=cg_k(:,1+(2*ib-2)*npw_k+ishf:(2*ib-1)*npw_k+ishf)
785 27271096 : cwavef(:,(ib-1)*npw_k+1:ib*npw_k,2)=cg_k(:,1+(2*ib-1)*npw_k+ishf:ib*2*npw_k+ishf)
786 : end do
787 : end if
788 : else
789 800 : ishf=(iblock-1)*npw_k*my_nspinor*blocksize
790 1600 : do ib=1,blocksize
791 : cwavef(:,(ib-1)*npw_k+1:ib*npw_k,ispinor_index)=&
792 392000 : & cg_k(:,1+(ib-1)*npw_k+ishf:ib*npw_k+ishf)
793 197200 : cwavef(:,(ib-1)*npw_k+1:ib*npw_k,jspinor_index)=zero
794 : end do
795 800 : call xmpi_sum(cwavef,mpi_enreg%comm_spinor,ierr)
796 : end if
797 : end if
798 :
799 79232 : if(ioption==1)then
800 : ! Multiplication by 2pi i (k+G)_alpha
801 7326 : gp2pi1=gprimd(alpha,1)*two_pi ; gp2pi2=gprimd(alpha,2)*two_pi ; gp2pi3=gprimd(alpha,3)*two_pi
802 7326 : kpt_cart=gp2pi1*dtset%kptns(1,ikpt)+gp2pi2*dtset%kptns(2,ikpt)+gp2pi3*dtset%kptns(3,ikpt)
803 1077648 : kg_k_cart_block(1:npw_k)=gp2pi1*kg_k(1,1:npw_k)+gp2pi2*kg_k(2,1:npw_k)+gp2pi3*kg_k(3,1:npw_k)+kpt_cart
804 7326 : if(gpu_cwavef) then
805 : #ifdef HAVE_OPENMP_OFFLOAD
806 : !$OMP TARGET TEAMS DISTRIBUTE MAP(to:kg_k_cart_block,cwavef) PRIVATE(ib)
807 : do ib=1,blocksize
808 : !$OMP PARALLEL DO PRIVATE(ipw,cwftmp)
809 : do ipw=1,npw_k
810 : cwftmp=-cwavef(2,ipw+(ib-1)*npw_k,1)*kg_k_cart_block(ipw)
811 : cwavef(2,ipw+(ib-1)*npw_k,1)=cwavef(1,ipw+(ib-1)*npw_k,1)*kg_k_cart_block(ipw)
812 : cwavef(1,ipw+(ib-1)*npw_k,1)=cwftmp
813 : end do
814 : end do
815 : if (my_nspinor==2) then
816 : !$OMP TARGET TEAMS DISTRIBUTE MAP(to:kg_k_cart_block,cwavef) PRIVATE(ib)
817 : do ib=1,blocksize
818 : !$OMP PARALLEL DO PRIVATE(ipw,cwftmp)
819 : do ipw=1,npw_k
820 : cwftmp=-cwavef(2,ipw+(ib-1)*npw_k,2)*kg_k_cart_block(ipw)
821 : cwavef(2,ipw+(ib-1)*npw_k,2)=cwavef(1,ipw+(ib-1)*npw_k,2)*kg_k_cart_block(ipw)
822 : cwavef(1,ipw+(ib-1)*npw_k,2)=cwftmp
823 : end do
824 : end do
825 : end if
826 : #endif
827 : else
828 39486 : do ib=1,blocksize
829 5730126 : do ipw=1,npw_k
830 5690640 : cwftmp=-cwavef(2,ipw+(ib-1)*npw_k,1)*kg_k_cart_block(ipw)
831 5690640 : cwavef(2,ipw+(ib-1)*npw_k,1)=cwavef(1,ipw+(ib-1)*npw_k,1)*kg_k_cart_block(ipw)
832 5690640 : cwavef(1,ipw+(ib-1)*npw_k,1)=cwftmp
833 5722800 : if (my_nspinor==2) then
834 3405600 : cwftmp=-cwavef(2,ipw+(ib-1)*npw_k,2)*kg_k_cart_block(ipw)
835 3405600 : cwavef(2,ipw+(ib-1)*npw_k,2)=cwavef(1,ipw+(ib-1)*npw_k,2)*kg_k_cart_block(ipw)
836 3405600 : cwavef(1,ipw+(ib-1)*npw_k,2)=cwftmp
837 : end if
838 : end do
839 : end do
840 : end if
841 71906 : else if(ioption==2)then
842 0 : ABI_ERROR("kinetic energy density tensor (taur_(alpha,beta)) is not yet implemented.")
843 : end if
844 :
845 : ! ---------- DMFT
846 79232 : if(paw_dmft%use_sc_dmft==1.or.paw_dmft%use_sc_dmft==10) then
847 : ! initialisation of DMFT arrays
848 3963256 : cwavef_rot(:,:,:,:) = zero
849 7260 : occ_diag(:) = zero
850 : ! occ_nd(:,:,:,:) = paw_dmft%occnd(:,:,:,ikpt,:)
851 :
852 : #ifdef HAVE_OPENMP_OFFLOAD
853 : !$OMP TARGET UPDATE FROM(cwavef) IF(gpu_cwavef)
854 : #endif
855 7260 : do ib=1,blocksize
856 3970076 : cwavef_rot(:, :, ib, :) = cwavef(:, 1+(ib-1)*npw_k:ib*npw_k, :)
857 : end do
858 : call rot_cg(paw_dmft%occnd(:,:,:,ikpt,isppol), cwavef_rot, npw_k, nband_k, blocksize,&
859 : & dtset%nspinor, paw_dmft%include_bands(1), paw_dmft%mbandc, occ_diag,&
860 220 : & (paw_dmft%dmft_solv == 6 .or. paw_dmft%dmft_solv == 7))
861 7260 : do ib=1,blocksize
862 3970076 : cwavef(:, 1+(ib-1)*npw_k:ib*npw_k, :) = cwavef_rot(:, :, ib, :)
863 : end do
864 : #ifdef HAVE_OPENMP_OFFLOAD
865 : !$OMP TARGET UPDATE TO(cwavef) IF(gpu_cwavef)
866 : #endif
867 :
868 7260 : occ_k(:) = occ_diag(:)
869 : end if
870 : ! ---------- END DMFT
871 :
872 79232 : call timab(538,1,tsec)
873 79232 : if (nspinor1TreatedByThisProc) then
874 : call prep_fourwf(rhoaug,blocksize,cwavef(:,:,1),wfraug,iblock,istwf_k,dtset%mgfft,mpi_enreg,&
875 : & nband_k,ndat,dtset%ngfft,npw_k,n4,n5,n6,occ_k,1,ucvol,&
876 78832 : & dtset%wtk(ikpt),l_nfft_blocks,gpu_option=gpu_option)
877 : end if
878 79232 : call timab(538,2,tsec)
879 95240 : if(dtset%nspinor==2)then
880 5656 : if (dtset%nspden==1) then
881 4856 : if (nspinor2TreatedByThisProc) then
882 : call prep_fourwf(rhoaug,blocksize,cwavef(:,:,2),wfraug,&
883 : & iblock,istwf_k,dtset%mgfft,mpi_enreg,&
884 : & nband_k,ndat,dtset%ngfft,npw_k,n4,n5,n6,occ_k,1,ucvol,&
885 4856 : & dtset%wtk(ikpt),l_nfft_blocks,gpu_option=gpu_option)
886 : end if
887 800 : else if(dtset%nspden==4 ) then
888 :
889 800 : if(gpu_option == ABI_GPU_KOKKOS) then
890 : #if defined HAVE_GPU && defined HAVE_YAKL
891 : ABI_MALLOC_MANAGED(cwavef_x,(/2,npw_k*blocksize/))
892 : ABI_MALLOC_MANAGED(cwavef_y,(/2,npw_k*blocksize/))
893 : #endif
894 : else
895 2400 : ABI_MALLOC(cwavef_x,(2,npw_k*blocksize))
896 1600 : ABI_MALLOC(cwavef_y,(2,npw_k*blocksize))
897 : end if
898 :
899 800 : if(gpu_cwavef) then
900 : #ifdef HAVE_OPENMP_OFFLOAD
901 : !$OMP TARGET ENTER DATA MAP(alloc:cwavef_x,cwavef_y)
902 : !$OMP TARGET TEAMS DISTRIBUTE MAP(to:cwavef,cwavef_x,cwavef_y) PRIVATE(ib)
903 : do ib=1,blocksize
904 : !$OMP PARALLEL DO PRIVATE(ipw)
905 : do ipw=1,npw_k
906 : cwavef_x(1,(ib-1)*npw_k+ipw)=cwavef(1,(ib-1)*npw_k+ipw,1)+cwavef(1,(ib-1)*npw_k+ipw,2)
907 : cwavef_x(2,(ib-1)*npw_k+ipw)=cwavef(2,(ib-1)*npw_k+ipw,1)+cwavef(2,(ib-1)*npw_k+ipw,2)
908 :
909 : cwavef_y(1,(ib-1)*npw_k+ipw)=cwavef(1,(ib-1)*npw_k+ipw,1)+cwavef(2,(ib-1)*npw_k+ipw,2)
910 : cwavef_y(2,(ib-1)*npw_k+ipw)=cwavef(2,(ib-1)*npw_k+ipw,1)-cwavef(1,(ib-1)*npw_k+ipw,2)
911 : end do
912 : end do
913 : #endif
914 : else
915 392000 : cwavef_x(:,:)=cwavef(:,:,1)+cwavef(:,:,2)
916 131200 : cwavef_y(1,:)=cwavef(1,:,1)+cwavef(2,:,2)
917 131200 : cwavef_y(2,:)=cwavef(2,:,1)-cwavef(1,:,2)
918 : end if
919 :
920 800 : call timab(538,1,tsec)
921 800 : if (nspinor1TreatedByThisProc) then
922 : call prep_fourwf(rhoaug_down,blocksize,cwavef(:,:,2),wfraug,&
923 : & iblock,istwf_k,dtset%mgfft,mpi_enreg,&
924 : & nband_k,ndat,dtset%ngfft,npw_k,n4,n5,n6,occ_k,1,ucvol,&
925 400 : & dtset%wtk(ikpt),l_nfft_blocks,gpu_option=gpu_option)
926 : end if
927 800 : if (nspinor2TreatedByThisProc) then
928 : call prep_fourwf(rhoaug_mx,blocksize,cwavef_x,wfraug,&
929 : & iblock,istwf_k,dtset%mgfft,mpi_enreg,&
930 : & nband_k,ndat,dtset%ngfft,npw_k,n4,n5,n6,occ_k,1,ucvol,&
931 400 : & dtset%wtk(ikpt),l_nfft_blocks,gpu_option=gpu_option)
932 : call prep_fourwf(rhoaug_my,blocksize,cwavef_y,wfraug,&
933 : & iblock,istwf_k,dtset%mgfft,mpi_enreg,&
934 : & nband_k,ndat,dtset%ngfft,npw_k,n4,n5,n6,occ_k,1,ucvol,&
935 400 : & dtset%wtk(ikpt),l_nfft_blocks,gpu_option=gpu_option)
936 : end if
937 800 : call timab(538,2,tsec)
938 :
939 800 : if(gpu_option == ABI_GPU_KOKKOS) then
940 : #if defined HAVE_GPU && defined HAVE_YAKL
941 : ABI_FREE_MANAGED(cwavef_x)
942 : ABI_FREE_MANAGED(cwavef_y)
943 : #endif
944 : else
945 : #ifdef HAVE_OPENMP_OFFLOAD
946 : !$OMP TARGET EXIT DATA MAP(delete:cwavef_x,cwavef_y) IF(gpu_cwavef)
947 : #endif
948 800 : ABI_FREE(cwavef_x)
949 800 : ABI_FREE(cwavef_y)
950 : end if
951 :
952 : end if
953 : end if
954 : end do !iblock
955 16008 : if(ioption==1) then
956 1206 : ABI_FREE(kg_k_cart_block)
957 : end if
958 :
959 16008 : if (associated(cwavef)) then
960 16008 : if(gpu_option == ABI_GPU_KOKKOS) then
961 : #if defined HAVE_GPU && defined HAVE_YAKL
962 : ABI_FREE_MANAGED(cwavef)
963 : #endif
964 : else
965 : #ifdef HAVE_OPENMP_OFFLOAD
966 : !$OMP TARGET EXIT DATA MAP(delete:cwavef) IF(gpu_cwavef)
967 : #endif
968 16008 : ABI_FREE(cwavef)
969 : end if
970 : end if
971 :
972 16008 : ABI_FREE(occ_k)
973 : #ifdef HAVE_OPENMP_OFFLOAD
974 : !$OMP TARGET EXIT DATA MAP(delete:cg_k) if(gpu_cwavef .and. transfer_cg)
975 : #endif
976 : end if ! paral_kgb
977 :
978 173694 : ABI_FREE(gbound)
979 :
980 173694 : if(gpu_option == ABI_GPU_KOKKOS) then
981 : #if defined HAVE_GPU && defined HAVE_YAKL
982 : ABI_FREE_MANAGED(kg_k)
983 : #endif
984 : else
985 173694 : ABI_FREE(kg_k)
986 : end if
987 :
988 173694 : bdtot_index=bdtot_index+nband_k
989 :
990 194675 : if (dtset%mkmem/=0) then
991 173694 : icg=icg+npw_k*my_nspinor*mband_mem !iband_me
992 173694 : ikg=ikg+npw_k
993 : end if
994 :
995 : end do ! ikpt
996 :
997 : #ifdef HAVE_OPENMP_OFFLOAD
998 : !$OMP TARGET UPDATE FROM(rhoaug) IF(gpu_option==ABI_GPU_OPENMP)
999 : #endif
1000 20981 : if(mpi_enreg%paral_kgb == 1) then
1001 4488 : call bandfft_kpt_set_ikpt(-1,mpi_enreg)
1002 4488 : if (dtset%nspden==4) then
1003 : ! Sum the contribution of the band and of the FFT
1004 4 : call xmpi_sum(rhoaug ,mpi_enreg%comm_bandspinorfft, ierr)
1005 4 : call xmpi_sum(rhoaug_down,mpi_enreg%comm_bandspinorfft, ierr)
1006 4 : call xmpi_sum(rhoaug_mx ,mpi_enreg%comm_bandspinorfft, ierr)
1007 4 : call xmpi_sum(rhoaug_my ,mpi_enreg%comm_bandspinorfft, ierr)
1008 39300 : rhoaug_up(:,:,:) = rhoaug(:,:,:)
1009 : else
1010 4484 : call xmpi_sum(rhoaug,mpi_enreg%comm_bandspinorfft,ierr)
1011 : end if
1012 : end if
1013 :
1014 : ! Transfer density on augmented fft grid to normal fft grid in real space
1015 : ! Take also into account the spin, to place it correctly in rhor.
1016 39779 : if(dtset%nspden==1 .or. dtset%nspden==2) then
1017 19024 : call fftpac(isppol,mpi_enreg,dtset%nspden,n1,n2,n3,n4,n5,n6,dtset%ngfft,rhor,rhoaug,1)
1018 1957 : else if(dtset%nspden==4) then
1019 1957 : ispden=1
1020 1957 : call fftpac(ispden,mpi_enreg,dtset%nspden,n1,n2,n3,n4,n5,n6,dtset%ngfft,rhor,rhoaug_up,1)
1021 1957 : ispden=2
1022 1957 : call fftpac(ispden,mpi_enreg,dtset%nspden,n1,n2,n3,n4,n5,n6,dtset%ngfft,rhor,rhoaug_mx,1)
1023 1957 : ispden=3
1024 1957 : call fftpac(ispden,mpi_enreg,dtset%nspden,n1,n2,n3,n4,n5,n6,dtset%ngfft,rhor,rhoaug_my,1)
1025 1957 : ispden=4
1026 1957 : call fftpac(ispden,mpi_enreg,dtset%nspden,n1,n2,n3,n4,n5,n6,dtset%ngfft,rhor,rhoaug_down,1)
1027 1957 : if(gpu_option == ABI_GPU_KOKKOS) then
1028 : #if defined HAVE_GPU && defined HAVE_YAKL
1029 : ABI_FREE_MANAGED(rhoaug_up)
1030 : ABI_FREE_MANAGED(rhoaug_down)
1031 : ABI_FREE_MANAGED(rhoaug_mx)
1032 : ABI_FREE_MANAGED(rhoaug_my)
1033 : #endif
1034 : else
1035 1957 : ABI_FREE(rhoaug_up)
1036 1957 : ABI_FREE(rhoaug_down)
1037 1957 : ABI_FREE(rhoaug_mx)
1038 1957 : ABI_FREE(rhoaug_my)
1039 : end if
1040 : end if
1041 :
1042 : end do ! isppol
1043 :
1044 18798 : if(gpu_option == ABI_GPU_KOKKOS) then
1045 : #if defined HAVE_GPU && defined HAVE_YAKL
1046 : if(associated(cwavef)) then
1047 : ABI_FREE_MANAGED(cwavef)
1048 : end if
1049 : if(associated(cwavefb)) then
1050 : ABI_FREE_MANAGED(cwavefb)
1051 : end if
1052 : ABI_FREE_MANAGED(rhoaug)
1053 : ABI_FREE_MANAGED(wfraug)
1054 : #endif
1055 : else
1056 18798 : if(associated(cwavef)) then
1057 14310 : ABI_FREE(cwavef)
1058 : end if
1059 18798 : if(associated(cwavefb)) then
1060 18798 : ABI_FREE(cwavefb)
1061 : endif
1062 : #ifdef HAVE_OPENMP_OFFLOAD
1063 : !$OMP TARGET EXIT DATA MAP(delete:rhoaug) IF(gpu_option==ABI_GPU_OPENMP)
1064 : #endif
1065 18798 : ABI_FREE(rhoaug)
1066 18798 : ABI_FREE(wfraug)
1067 : end if
1068 :
1069 18798 : if(allocated(cwavef_rot)) then
1070 44 : ABI_FREE(cwavef_rot)
1071 44 : ABI_FREE(occ_diag)
1072 : ! ABI_FREE(occ_nd)
1073 : end if
1074 :
1075 : ! Recreate full rhor on all proc.
1076 18798 : call timab(48,1,tsec)
1077 18798 : call timab(71,1,tsec)
1078 18798 : spaceComm=mpi_enreg%comm_cell
1079 18798 : if (mpi_enreg%paral_hf==1)spaceComm=mpi_enreg%comm_kpt
1080 18798 : if(mpi_enreg%paral_kgb==1)spaceComm=mpi_enreg%comm_kpt
1081 18798 : call xmpi_sum(rhor,spaceComm,ierr)
1082 18798 : call timab(71,2,tsec)
1083 18798 : call timab(48,2,tsec)
1084 :
1085 18798 : call timab(799,2,tsec)
1086 18798 : call timab(549,1,tsec)
1087 :
1088 37596 : if(ioption==1 .or. ioption==2) then
1089 : !$OMP PARALLEL DO COLLAPSE(2)
1090 5664 : do ispden=1,dtset%nspden
1091 116616666 : do ifft=1,dtset%nfft
1092 116614194 : taur_alphabeta(ifft,ispden,alpha,beta) = rhor(ifft,ispden)
1093 : end do
1094 : end do
1095 : end if
1096 :
1097 : end do ! beta=1,nbeta
1098 : end do ! alpha=1,nalpha
1099 :
1100 : !Compute the trace over the kinetic energy density tensor. i.e. Sum of the 3 diagonal elements.
1101 17150 : if(ioption==1)then
1102 : ! zero rhor array in real space
1103 1888 : do ispden=1,dtset%nspden
1104 : !$OMP PARALLEL DO
1105 38872222 : do ifft=1,dtset%nfft
1106 38871398 : rhor(ifft,ispden)=zero
1107 : end do
1108 : end do
1109 3296 : do alpha = 1, nalpha
1110 : !$OMP PARALLEL DO COLLAPSE(2)
1111 6488 : do ispden=1,dtset%nspden
1112 116616666 : do ifft=1,dtset%nfft
1113 116614194 : rhor(ifft,ispden) = rhor(ifft,ispden) + taur_alphabeta(ifft,ispden,alpha,1)
1114 : end do
1115 : end do
1116 : end do
1117 : end if
1118 :
1119 17150 : nfftot=dtset%ngfft(1) * dtset%ngfft(2) * dtset%ngfft(3)
1120 :
1121 : !Add extfpmd electrons contributions to density on coarse grid.
1122 : !When using a fine grid, space-dependant contributions to the
1123 : !density are added in the pawmkrho subroutine.
1124 17150 : if(present(extfpmd)) then
1125 8295 : if(associated(extfpmd)) then
1126 0 : if(extfpmd%version==10.and.allocated(extfpmd%nelectarr)) then
1127 0 : rhor(:,:)=rhor(:,:)+extfpmd%nelectarr(:,:)/ucvol/dtset%nspden
1128 : else
1129 0 : rhor(:,:)=rhor(:,:)+extfpmd%nelect/ucvol/dtset%nspden
1130 : end if
1131 : end if
1132 : end if
1133 :
1134 17150 : select case (ioption)
1135 : case (0, 1)
1136 : call symrhg(1,gprimd,irrzon,mpi_enreg,dtset%nfft,nfftot,dtset%ngfft,dtset%nspden,dtset%nsppol,dtset%nsym,&
1137 17150 : phnons,rhog,rhor,rprimd,dtset%symafm,dtset%symrel,dtset%tnons)
1138 17150 : if(ioption==1)then
1139 : !$OMP PARALLEL DO
1140 31356088 : do ifft=1,dtset%nfft
1141 70225598 : do ispden=1,dtset%nspden
1142 70225598 : rhor(ifft,ispden)=1.0d0/2.0d0*rhor(ifft,ispden)
1143 : end do
1144 94066616 : rhog(:,ifft)=1.0d0/2.0d0*rhog(:,ifft)
1145 : end do
1146 : end if
1147 : case (2)
1148 17150 : ABI_BUG('kinetic energy density tensor (taur_(alpha,beta)) is not yet implemented.')
1149 : !call symtaug(1,gprimd,irrzon,mpi_enreg,dtset%nfft,nfftot,dtset%ngfft,dtset%nspden,dtset%nsppol,dtset%nsym,&
1150 : !dtset%paral_kgb,phnons,rhog,rhor,rprimd,dtset%symafm,dtset%symrel)
1151 : end select
1152 :
1153 17150 : call timab(549,2,tsec)
1154 :
1155 : !We now have both rho(r) and rho(G), symmetrized, and if dtset%nsppol=2
1156 : !we also have the spin-up density, symmetrized, in rhor(:,2).
1157 : !In case of non collinear magnetism, we have rho,mx,my,mz. No symmetry is applied
1158 :
1159 17150 : call timab(799,1,tsec)
1160 :
1161 17150 : if(ioption==1 .or. ioption==2) then
1162 824 : ABI_FREE(taur_alphabeta)
1163 : end if
1164 :
1165 : !Find and print minimum and maximum total electron density
1166 : !(or total kinetic energy density, or total element of kinetic energy density tensor) and locations
1167 17150 : if (l_printout) then
1168 15945 : call wrtout(std_out,' mkrho: echo density (plane-wave part only)','COLL')
1169 15945 : call prtrhomxmn(std_out,mpi_enreg,dtset%nfft,dtset%ngfft,dtset%nspden,1,rhor,optrhor=ioption,ucvol=ucvol)
1170 : end if
1171 :
1172 17150 : call timab(799,2,tsec)
1173 17150 : call timab(790+tim_mkrho,2,tsec)
1174 :
1175 : DBG_EXIT("COLL")
1176 :
1177 17150 : end subroutine mkrho
1178 : !!***
1179 :
1180 : !!****f* m_mkrho/initro
1181 : !! NAME
1182 : !! initro
1183 : !!
1184 : !! FUNCTION
1185 : !! Initialize the density using either:
1186 : !! - a gaussian of adjustable decay length (norm-conserving psp)
1187 : !! - PS atomic valence density from psp file (PAW or NC psps with valence change in the pp file)
1188 : !!
1189 : !! INPUTS
1190 : !! atindx(natom)=index table for atoms (see gstate.f)
1191 : !! densty(ntypat,4)=parameters for initialisation of the density of each atom type
1192 : !! gmet(3,3)=reciprocal space metric (Bohr**-2)
1193 : !! gsqcut=cutoff G**2 for included G s in fft box (larger sphere).
1194 : !! izero=if 1, unbalanced components of rho(g) have to be set to zero
1195 : !! mgfft=maximum size of 1D FFTs
1196 : !! mpi_enreg=information about mpi parallelization
1197 : !! mqgrid=number of grid pts in q array for n^AT(q) spline.
1198 : !! natom=number of atoms in cell.
1199 : !! nattyp(ntypat)=number of atoms of each type in cell.
1200 : !! nfft=(effective) number of FFT grid points (for this processor)
1201 : !! ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft
1202 : !! ntypat=number of types of atoms in cell.
1203 : !! nspden=number of spin-density components
1204 : !! psps<type(pseudopotential_type)>=variables related to pseudopotentials
1205 : !! pawtab(ntypat*usepaw) <type(pawtab_type)>=paw tabulated starting data
1206 : !! ph1d(2,3*(2*mgfft+1)*natom)=1-dim phase information for given atom coordinates.
1207 : !! qgrid(mqgrid)=q grid for spline atomic valence density n^AT(q) from 0 to qmax.
1208 : !! spinat(3,natom)=initial spin of each atom, in unit of hbar/2.
1209 : !! ucvol=unit cell volume (Bohr**3).
1210 : !! usepaw= 0 for non paw calculation; =1 for paw calculation
1211 : !! zion(ntypat)=charge on each type of atom (real number)
1212 : !! znucl(ntypat)=atomic number, for each type of atom
1213 : !!
1214 : !! OUTPUT
1215 : !! rhog(2,nfft)=initialized total density in reciprocal space
1216 : !! rhor(nfft,nspden)=initialized total density in real space.
1217 : !! as well as spin-up part if spin-polarized
1218 : !!
1219 : !! SOURCE
1220 :
1221 3067 : subroutine initro(atindx,densty,gmet,gsqcut,izero,mgfft,mpi_enreg,mqgrid,natom,nattyp,&
1222 3067 : & nfft,ngfft,nspden,ntypat,psps,pawtab,ph1d,qgrid,rhog,rhor,spinat,ucvol,usepaw,zion,znucl)
1223 :
1224 : !Arguments ------------------------------------
1225 : !scalars
1226 : integer,intent(in) :: izero,mgfft,mqgrid,natom,nfft,nspden,ntypat
1227 : integer,intent(in) :: usepaw
1228 : real(dp),intent(in) :: gsqcut,ucvol
1229 : type(mpi_type),intent(in) :: mpi_enreg
1230 : type(pseudopotential_type),intent(in) :: psps
1231 : !arrays
1232 : integer,intent(in) :: atindx(natom),nattyp(ntypat),ngfft(18)
1233 : real(dp),intent(in) :: densty(ntypat,4),gmet(3,3),ph1d(2,3*(2*mgfft+1)*natom)
1234 : real(dp),intent(in) :: qgrid(mqgrid),spinat(3,natom),zion(ntypat)
1235 : real(dp),intent(in) :: znucl(ntypat)
1236 : real(dp),intent(out) :: rhog(2,nfft),rhor(nfft,nspden)
1237 : type(pawtab_type),intent(in) :: pawtab(ntypat*usepaw)
1238 :
1239 : !Local variables-------------------------------
1240 : !scalars
1241 : integer,parameter :: im=2,re=1
1242 : integer :: i1,i2,i3,ia,ia1,ia2,id1,id2,id3,ig1,ig2,ig3,ii,ispden
1243 : integer :: itypat,jj,jtemp,me_fft,n1,n2,n3,nproc_fft
1244 : real(dp),parameter :: tolfix=1.000000001_dp
1245 : real(dp) :: aa,alf2pi2,bb,cc,cutoff,dd,diff,dq,dq2div6,dqm1,fact,fact0,gmag
1246 : real(dp) :: gsquar,rhoat,sfi,sfr
1247 : real(dp) :: xnorm
1248 : character(len=500) :: message
1249 : !arrays
1250 3067 : integer, contiguous, pointer :: fftn2_distrib(:),ffti2_local(:),fftn3_distrib(:),ffti3_local(:)
1251 6134 : real(dp) :: length(ntypat)
1252 3067 : real(dp),allocatable :: work(:), spinat_indx(:,:)
1253 6134 : logical :: use_gaussian(ntypat)
1254 : ! *************************************************************************
1255 :
1256 3067 : if (nspden==4) then
1257 120 : ABI_COMMENT('initro: might work yet for nspden=4 (not checked)')
1258 : !write(std_out,*)' spinat',spinat(1:3,1:natom)
1259 : end if
1260 :
1261 3067 : n1=ngfft(1); n2=ngfft(2); n3=ngfft(3); me_fft=ngfft(11); nproc_fft=ngfft(10)
1262 :
1263 9201 : ABI_MALLOC(work,(nfft))
1264 9201 : ABI_MALLOC(spinat_indx,(3,natom))
1265 :
1266 : ! Get the distrib associated with this fft_grid
1267 3067 : call ptabs_fourdp(mpi_enreg,n2,n3,fftn2_distrib,ffti2_local,fftn3_distrib,ffti3_local)
1268 :
1269 : ! Transfer the spinat array to an array in which the atoms have the proper order, type by type.
1270 10769 : do ia=1,natom
1271 33875 : spinat_indx(:,atindx(ia))=spinat(:,ia)
1272 : end do
1273 :
1274 : ! Check whether the values of spinat are acceptable
1275 3067 : if (nspden==2)then
1276 385 : ia1=1
1277 828 : do itypat=1,ntypat
1278 : !ia1,ia2 sets range of loop over atoms:
1279 443 : ia2=ia1+nattyp(itypat)-1
1280 1082 : do ia=ia1,ia2
1281 639 : if( sqrt(spinat_indx(1,ia)**2+spinat_indx(2,ia)**2+spinat_indx(3,ia)**2) &
1282 443 : > abs(zion(itypat))*(1.0_dp + epsilon(0.0_dp)) ) then
1283 : write(message, '(a,i0,a,a,3es11.4,a,a,a,es11.4)' )&
1284 0 : ' For type-ordered atom number ',ia,ch10,&
1285 0 : ' input spinat=',spinat_indx(:,ia),' is larger, in magnitude,',ch10,&
1286 0 : ' than zion(ia)=',zion(itypat)
1287 0 : call wrtout([std_out, ab_out], message)
1288 : end if
1289 : end do
1290 828 : ia1=ia2+1
1291 : end do
1292 : end if
1293 :
1294 : ! Compute the decay length of each type of atom depending on data available in pseudos.
1295 3067 : jtemp=0
1296 7402 : do itypat=1,ntypat
1297 4335 : use_gaussian(itypat)=.true.
1298 4335 : if (usepaw==0) use_gaussian(itypat) = .not. psps%nctab(itypat)%has_tvale
1299 4335 : if (usepaw==1) use_gaussian(itypat)=(pawtab(itypat)%has_tvale==0)
1300 4335 : if (.not.use_gaussian(itypat)) jtemp=jtemp+1
1301 :
1302 7402 : if (use_gaussian(itypat)) then
1303 2916 : length(itypat) = atom_length(densty(itypat,1),zion(itypat),znucl(itypat))
1304 : write(message,'(a,i3,a,f12.4,a,a,a,f12.4,a,i3,a,es12.4,a)' )&
1305 2916 : ' initro: for itypat=',itypat,', take decay length=',length(itypat),',',ch10,&
1306 5832 : ' initro: indeed, coreel=',znucl(itypat)-zion(itypat),', nval=',int(zion(itypat)),' and densty=',densty(itypat,1),'.'
1307 2916 : call wrtout(std_out,message)
1308 : else
1309 1419 : write(message,"(a,i3,a)")' initro: for itypat=',itypat,", take pseudo charge density from pp file"
1310 1419 : call wrtout(std_out,message)
1311 : end if
1312 : end do
1313 :
1314 3067 : if (jtemp>0) then
1315 913 : dq=(qgrid(mqgrid)-qgrid(1))/dble(mqgrid-1)
1316 913 : dqm1=1.0_dp/dq
1317 913 : dq2div6=dq**2/6.0_dp
1318 : end if
1319 :
1320 3067 : cutoff=gsqcut*tolfix
1321 3067 : xnorm=1.0_dp/ucvol
1322 :
1323 3067 : id1=n1/2+2
1324 3067 : id2=n2/2+2
1325 3067 : id3=n3/2+2
1326 :
1327 3067 : if(nspden /= 4) then
1328 :
1329 6279 : do ispden=nspden,1,-1
1330 : ! This loop overs spins will actually be as follows:
1331 : ! ispden=2 for spin up
1332 : ! ispden=1 for total spin (also valid for non-spin-polarized calculations)
1333 : !
1334 : ! The reverse ispden order is chosen, in order to end up with
1335 : ! rhog containing the proper total density.
1336 287175098 : rhog = zero
1337 :
1338 3332 : ia1=1
1339 7958 : do itypat=1,ntypat
1340 4626 : if (use_gaussian(itypat)) alf2pi2=(two_pi*length(itypat))**2
1341 :
1342 : ! ia1,ia2 sets range of loop over atoms:
1343 4626 : ia2=ia1+nattyp(itypat)-1
1344 :
1345 : !$OMP PARALLEL DO &
1346 : !$OMP& PRIVATE(i3,i2,i1,ig3,ig2,ig1,ii,jj,gsquar,fact,sfr,sfi,fact0,rhoat,aa,bb,cc,dd,gmag,diff)
1347 116655 : do i3=1,n3
1348 112029 : ig3=i3-(i3/id3)*n3-1
1349 3378276 : do i2=1,n2
1350 3261621 : ig2=i2-(i2/id2)*n2-1
1351 3373650 : if (fftn2_distrib(i2)==me_fft) then
1352 140966977 : do i1=1,n1
1353 :
1354 137826164 : ig1=i1-(i1/id1)*n1-1
1355 137826164 : ii=i1+n1*(ffti2_local(i2)-1+(n2/nproc_fft)*(i3-1))
1356 : gsquar=dble(ig1*ig1)*gmet(1,1)+dble(ig2*ig2)*gmet(2,2)+&
1357 : dble(ig3*ig3)*gmet(3,3)+dble(2*ig1*ig2)*gmet(1,2)+&
1358 137826164 : dble(2*ig2*ig3)*gmet(2,3)+dble(2*ig3*ig1)*gmet(3,1)
1359 :
1360 : ! Skip G**2 outside cutoff:
1361 140966977 : if (gsquar<=cutoff) then
1362 :
1363 : ! Assemble structure factor over all atoms of given type,
1364 : ! also taking into account the spin-charge on each atom:
1365 :
1366 56338141 : sfr=zero;sfi=zero
1367 56338141 : if (ispden==1) then
1368 149015914 : do ia=ia1,ia2
1369 103544853 : sfr=sfr+phre_ini(ig1,ig2,ig3,ia)
1370 149015914 : sfi=sfi-phimag_ini(ig1,ig2,ig3,ia)
1371 : end do
1372 45471061 : if (use_gaussian(itypat)) then
1373 25154892 : sfr=sfr*zion(itypat)
1374 25154892 : sfi=sfi*zion(itypat)
1375 : end if
1376 : else
1377 10867080 : fact0=half;if (.not.use_gaussian(itypat)) fact0=half/zion(itypat)
1378 21697592 : do ia=ia1,ia2
1379 : ! Here, take care only of the z component
1380 10830512 : fact=fact0*(zion(itypat)+spinat_indx(3,ia))
1381 10830512 : sfr=sfr+phre_ini(ig1,ig2,ig3,ia)*fact
1382 21697592 : sfi=sfi-phimag_ini(ig1,ig2,ig3,ia)*fact
1383 : end do
1384 : end if
1385 :
1386 : ! Charge density integrating to one
1387 56338141 : if (use_gaussian(itypat)) then
1388 34433254 : rhoat=xnorm*exp(-gsquar*alf2pi2)
1389 : ! Multiply structure factor times rhoat (atomic density in reciprocal space)
1390 34433254 : rhog(re,ii)=rhog(re,ii)+sfr*rhoat
1391 34433254 : rhog(im,ii)=rhog(im,ii)+sfi*rhoat
1392 : else
1393 21904887 : gmag=sqrt(gsquar)
1394 21904887 : jj=1+int(gmag*dqm1)
1395 21904887 : diff=gmag-qgrid(jj)
1396 21904887 : bb = diff*dqm1
1397 21904887 : aa = one-bb
1398 21904887 : cc = aa*(aa**2-one)*dq2div6
1399 21904887 : dd = bb*(bb**2-one)*dq2div6
1400 21904887 : if (usepaw == 1) then
1401 : rhoat=(aa*pawtab(itypat)%tvalespl(jj,1)+bb*pawtab(itypat)%tvalespl(jj+1,1)+&
1402 16650313 : & cc*pawtab(itypat)%tvalespl(jj,2)+dd*pawtab(itypat)%tvalespl(jj+1,2)) *xnorm
1403 5254574 : else if (usepaw == 0) then
1404 : rhoat=(aa*psps%nctab(itypat)%tvalespl(jj,1)+bb*psps%nctab(itypat)%tvalespl(jj+1,1)+&
1405 5254574 : cc*psps%nctab(itypat)%tvalespl(jj,2)+dd*psps%nctab(itypat)%tvalespl(jj+1,2))*xnorm
1406 : else
1407 0 : ABI_BUG('Initialization of density is non consistent.')
1408 : end if
1409 : ! Multiply structure factor times rhoat (atomic density in reciprocal space)
1410 21904887 : rhog(re,ii)=rhog(re,ii)+sfr*rhoat
1411 21904887 : rhog(im,ii)=rhog(im,ii)+sfi*rhoat
1412 : end if
1413 :
1414 : end if
1415 :
1416 : end do ! i1
1417 : end if
1418 : end do ! i2
1419 : end do ! i3
1420 7958 : ia1=ia2+1
1421 :
1422 : end do ! itypat
1423 :
1424 : ! Set contribution of unbalanced components to zero
1425 3332 : if (izero==1) then
1426 951 : call zerosym(rhog,2,n1,n2,n3,comm_fft=mpi_enreg%comm_fft,distribfft=mpi_enreg%distribfft)
1427 : end if
1428 : !write(std_out,*)"initro: ispden, ucvol * rhog(:2,1)",ispden, ucvol * rhog(:2,1)
1429 :
1430 : ! Note, we end with ispden=1, so that rhog contains the total density
1431 3332 : call fourdp(1,rhog,work,1,mpi_enreg,nfft,1,ngfft,0)
1432 95730201 : rhor(:,ispden)=work(:)
1433 : end do ! End loop on spins
1434 :
1435 : else if (nspden==4) then
1436 :
1437 600 : do ispden=nspden,1,-1
1438 : ! This loop overs spins will actually be as follows:
1439 : ! ispden=2,3,4 for mx,my,mz
1440 : ! ispden=1 for total spin (also valid for non-spin-polarized calculations)
1441 : ! The reverse ispden order is chosen, in order to end up with
1442 : ! rhog containing the proper total density.
1443 :
1444 18974448 : rhog(:,:)=zero
1445 :
1446 480 : ia1=1
1447 1088 : do itypat=1,ntypat
1448 :
1449 608 : if (use_gaussian(itypat)) alf2pi2=(two_pi*length(itypat))**2
1450 :
1451 : ! ia1,ia2 sets range of loop over atoms:
1452 608 : ia2=ia1+nattyp(itypat)-1
1453 608 : ii=0
1454 14612 : do i3=1,n3
1455 14004 : ig3=i3-(i3/id3)*n3-1
1456 335360 : do i2=1,n2
1457 320748 : ig2=i2-(i2/id2)*n2-1
1458 334752 : if (fftn2_distrib(i2)==me_fft) then
1459 8449512 : do i1=1,n1
1460 :
1461 8137980 : ig1=i1-(i1/id1)*n1-1
1462 8137980 : ii=ii+1
1463 : gsquar=dble(ig1*ig1)*gmet(1,1)+dble(ig2*ig2)*gmet(2,2)+&
1464 : dble(ig3*ig3)*gmet(3,3)+dble(2*ig1*ig2)*gmet(1,2)+&
1465 8137980 : dble(2*ig2*ig3)*gmet(2,3)+dble(2*ig3*ig1)*gmet(3,1)
1466 :
1467 : ! Skip G**2 outside cutoff:
1468 8449512 : if (gsquar<=cutoff) then
1469 :
1470 : ! Assemble structure factor over all atoms of given type,
1471 : ! also taking into account the spin-charge on each atom:
1472 2750480 : sfr=zero;sfi=zero
1473 2750480 : if(ispden==1)then
1474 1989487 : do ia=ia1,ia2
1475 1301867 : sfr=sfr+phre_ini(ig1,ig2,ig3,ia)
1476 1989487 : sfi=sfi-phimag_ini(ig1,ig2,ig3,ia)
1477 : end do
1478 687620 : if (use_gaussian(itypat)) then
1479 420048 : sfr=sfr*zion(itypat)
1480 420048 : sfi=sfi*zion(itypat)
1481 : end if
1482 : else
1483 2062860 : fact0=one;if (.not.use_gaussian(itypat)) fact0=one/zion(itypat)
1484 5968461 : do ia=ia1,ia2
1485 : ! Here, take care of the components of m
1486 3905601 : fact=fact0*spinat_indx(ispden-1,ia)
1487 3905601 : sfr=sfr+phre_ini(ig1,ig2,ig3,ia)*fact
1488 5968461 : sfi=sfi-phimag_ini(ig1,ig2,ig3,ia)*fact
1489 : end do
1490 : end if
1491 :
1492 : ! Charge density integrating to one
1493 2750480 : if (use_gaussian(itypat)) then
1494 1680192 : rhoat=xnorm*exp(-gsquar*alf2pi2)
1495 : else
1496 1070288 : gmag=sqrt(gsquar)
1497 1070288 : jj=1+int(gmag*dqm1)
1498 1070288 : diff=gmag-qgrid(jj)
1499 1070288 : bb = diff*dqm1
1500 1070288 : aa = one-bb
1501 1070288 : cc = aa*(aa**2-one)*dq2div6
1502 1070288 : dd = bb*(bb**2-one)*dq2div6
1503 1070288 : if (usepaw == 1) then
1504 : rhoat=(aa*pawtab(itypat)%tvalespl(jj,1)+bb*pawtab(itypat)%tvalespl(jj+1,1)+&
1505 615284 : & cc*pawtab(itypat)%tvalespl(jj,2)+dd*pawtab(itypat)%tvalespl(jj+1,2)) *xnorm
1506 455004 : else if (usepaw == 0) then
1507 : rhoat=(aa*psps%nctab(itypat)%tvalespl(jj,1)+bb*psps%nctab(itypat)%tvalespl(jj+1,1)+&
1508 455004 : cc*psps%nctab(itypat)%tvalespl(jj,2)+dd*psps%nctab(itypat)%tvalespl(jj+1,2))*xnorm
1509 : else
1510 0 : ABI_BUG('Initialization of density is non consistent.')
1511 : end if
1512 : end if
1513 :
1514 : ! Multiply structure factor times rhoat (atomic density in reciprocal space)
1515 2750480 : rhog(re,ii)=rhog(re,ii)+sfr*rhoat
1516 2750480 : rhog(im,ii)=rhog(im,ii)+sfi*rhoat
1517 : end if
1518 :
1519 : end do ! i1
1520 : end if
1521 : end do ! i2
1522 : end do ! i3
1523 1088 : ia1=ia2+1
1524 : end do ! itypat
1525 :
1526 : ! Set contribution of unbalanced components to zero
1527 480 : if (izero==1) then
1528 180 : call zerosym(rhog,2,n1,n2,n3,comm_fft=mpi_enreg%comm_fft,distribfft=mpi_enreg%distribfft)
1529 : end if
1530 : !write(std_out,*)"initro: ispden, ucvol * rhog(:2,1)",ispden, ucvol * rhog(:2,1)
1531 :
1532 : ! Note, we end with ispden=1, so that rhog contains the total density
1533 480 : call fourdp(1,rhog,work,1,mpi_enreg,nfft,1,ngfft,0)
1534 6325256 : rhor(:,ispden)=work(:)
1535 : end do ! ispden
1536 :
1537 : ! Non-collinear magnetism: avoid zero magnetization, because it produces numerical instabilities
1538 : ! Add a small real to the magnetization
1539 54917 : if (all(abs(spinat(:,:))<tol10)) rhor(:,4)=rhor(:,4)+tol14
1540 :
1541 : end if ! nspden==4
1542 :
1543 3067 : ABI_FREE(spinat_indx)
1544 6134 : ABI_FREE(work)
1545 :
1546 : contains
1547 :
1548 : !Real and imaginary parts of phase.
1549 119582833 : function phr_ini(x1,y1,x2,y2,x3,y3)
1550 :
1551 : real(dp) :: phr_ini
1552 : real(dp),intent(in) :: x1,x2,x3,y1,y2,y3
1553 119582833 : phr_ini=(x1*x2-y1*y2)*x3-(y1*x2+x1*y2)*y3
1554 : end function phr_ini
1555 :
1556 119582833 : function phi_ini(x1,y1,x2,y2,x3,y3)
1557 :
1558 : real(dp) :: phi_ini
1559 : real(dp),intent(in) :: x1,x2,x3,y1,y2,y3
1560 119582833 : phi_ini=(x1*x2-y1*y2)*y3+(y1*x2+x1*y2)*x3
1561 : end function phi_ini
1562 :
1563 478331332 : function ph1_ini(nri,ig1,ia)
1564 :
1565 : real(dp) :: ph1_ini
1566 : integer,intent(in) :: nri,ig1,ia
1567 478331332 : ph1_ini=ph1d(nri,ig1+1+n1+(ia-1)*(2*n1+1))
1568 478331332 : end function ph1_ini
1569 :
1570 478331332 : function ph2_ini(nri,ig2,ia)
1571 :
1572 : real(dp) :: ph2_ini
1573 : integer,intent(in) :: nri,ig2,ia
1574 478331332 : ph2_ini=ph1d(nri,ig2+1+n2+(ia-1)*(2*n2+1)+natom*(2*n1+1))
1575 478331332 : end function ph2_ini
1576 :
1577 478331332 : function ph3_ini(nri,ig3,ia)
1578 :
1579 : real(dp) :: ph3_ini
1580 : integer,intent(in) :: nri,ig3,ia
1581 478331332 : ph3_ini=ph1d(nri,ig3+1+n3+(ia-1)*(2*n3+1)+natom*(2*n1+1+2*n2+1))
1582 478331332 : end function ph3_ini
1583 :
1584 119582833 : function phre_ini(ig1,ig2,ig3,ia)
1585 :
1586 : real(dp) :: phre_ini
1587 : integer,intent(in) :: ig1,ig2,ig3,ia
1588 : phre_ini=phr_ini(ph1_ini(re,ig1,ia),ph1_ini(im,ig1,ia),&
1589 119582833 : & ph2_ini(re,ig2,ia),ph2_ini(im,ig2,ia),ph3_ini(re,ig3,ia),ph3_ini(im,ig3,ia))
1590 119582833 : end function phre_ini
1591 :
1592 119582833 : function phimag_ini(ig1,ig2,ig3,ia)
1593 :
1594 : real(dp) :: phimag_ini
1595 : integer,intent(in) :: ig1,ig2,ig3,ia
1596 : phimag_ini=phi_ini(ph1_ini(re,ig1,ia),ph1_ini(im,ig1,ia),&
1597 119582833 : & ph2_ini(re,ig2,ia),ph2_ini(im,ig2,ia),ph3_ini(re,ig3,ia),ph3_ini(im,ig3,ia))
1598 119582833 : end function phimag_ini
1599 :
1600 : end subroutine initro
1601 : !!***
1602 :
1603 : !!****f* m_mkrho/prtrhomxmn
1604 : !! NAME
1605 : !! prtrhomxmn
1606 : !!
1607 : !! FUNCTION
1608 : !! If option==1, compute the maximum and minimum of the density (and spin-polarization if nspden==2), and print it.
1609 : !! If option==2, also compute and print the second maximum or minimum
1610 : !!
1611 : !! INPUTS
1612 : !! iout=unit for output file
1613 : !! mpi_enreg=information about MPI parallelization
1614 : !! nfft=(effective) number of FFT grid points (for this processor)
1615 : !! ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft
1616 : !! nspden=number of spin-density components
1617 : !! option, see above
1618 : !! optrhor=option for rhor (If optrhor==0, rhor is expected to be the electron density)
1619 : !! (If optrhor==1, rhor is expected to be the kinetic energy density (taur))
1620 : !! (If optrhor==2, rhor is expected to be the gradient of the electron density (grhor))
1621 : !! (If optrhor==3, rhor is expected to be the laplacian of the electron density (lrhor))
1622 : !! (If optrhor==4, rhor is expected to be the ELF (elfr))
1623 : !! rhor(nfft,nspden)=electron density (electrons/bohr^3)
1624 : !!
1625 : !! NOTES
1626 : !! The tolerance tol12 aims at giving a machine-independent ordering.
1627 : !! (this trick is used in bonds.f, listkk.f, prtrhomxmn.f and rsiaf9.f)
1628 : !!
1629 : !! SOURCE
1630 :
1631 68932 : subroutine prtrhomxmn(iout,mpi_enreg,nfft,ngfft,nspden,option,rhor,optrhor,ucvol)
1632 :
1633 : !Arguments ------------------------------------
1634 : !scalars
1635 : integer,intent(in) :: iout,nfft,nspden,option
1636 : type(MPI_type),intent(in) :: mpi_enreg
1637 : integer,intent(in),optional :: optrhor
1638 : real(dp),intent(in),optional :: ucvol
1639 : !arrays
1640 : integer,intent(in) :: ngfft(18)
1641 : real(dp),intent(in) :: rhor(nfft,nspden)
1642 :
1643 : !Local variables-------------------------------
1644 : !scalars
1645 : integer :: i1,i2,i3,ierr,ifft,ii,iisign,iitems,index1,ioptrhor
1646 : integer :: index2,indsign,iproc,istart,me,n1,n2,n3,nitems
1647 : integer :: nfft_,nfftot,nproc,spaceComm
1648 : real(dp) :: temp,value1,value2
1649 : character(len=500) :: message,txt1_in_mssg,txt2_in_mssg,txt3_in_mssg
1650 : logical :: reduce=.false.
1651 : !arrays
1652 68932 : integer,allocatable :: iindex(:,:,:),index_fft(:,:,:,:)
1653 : real(dp) :: rhomn1(4),rhomn2(4),rhomx1(4),rhomx2(4),ri_rhomn1(3,4)
1654 : real(dp) :: ri_rhomn2(3,4),ri_rhomx1(3,4),ri_rhomx2(3,4),ri_zetmn1(3,2)
1655 : real(dp) :: ri_zetmn2(3,2),ri_zetmx1(3,2),ri_zetmx2(3,2),zetmn1(2)
1656 : real(dp) :: zetmn2(2),zetmx1(2),zetmx2(2)
1657 68932 : real(dp),allocatable :: array(:),coord(:,:,:,:),value(:,:,:),integrated(:)
1658 68932 : real(dp),allocatable :: value_fft(:,:,:)
1659 : ! *************************************************************************
1660 :
1661 68932 : if(.not.(present(optrhor))) then
1662 : ioptrhor=0
1663 : else
1664 16005 : ioptrhor=optrhor
1665 : end if
1666 :
1667 68932 : if(option/=1 .and. option/=2)then
1668 0 : ABI_BUG(sjoin(' Option must be 1 or 2, while it is:', itoa(option)))
1669 : end if
1670 :
1671 68932 : if (mpi_enreg%nproc_wvl>1) then
1672 : ! nfft is always the potential size (in GGA, the density has buffers).
1673 0 : nfft_ = ngfft(1) * ngfft(2) * mpi_enreg%nscatterarr(mpi_enreg%me_wvl, 2)
1674 0 : n1 = ngfft(1)
1675 0 : n2 = ngfft(2)
1676 0 : n3 = sum(mpi_enreg%nscatterarr(:, 2))
1677 0 : istart = mpi_enreg%nscatterarr(mpi_enreg%me_wvl, 4)
1678 : else
1679 68932 : nfft_ = nfft
1680 68932 : n1=ngfft(1) ; n2=ngfft(2) ; n3=ngfft(3)
1681 68932 : istart = 0
1682 : end if
1683 :
1684 : !--------------------------------------------------------------------------
1685 : !One has to determine the maximum and minimum (etc...) values
1686 : !over all space, and then output it, as well as to identify
1687 : !the point at which it occurs ...
1688 : !This will require a bit of data exchange, and correct indirect indexing ...
1689 :
1690 : !For the local processor, find different items :
1691 : !maximum and minimum total electron density and locations
1692 : !and also spin-polarisation and magnetization
1693 : !also keep the second maximal or minimal value
1694 : if(nspden==1)nitems=1 ! Simply the total density
1695 68932 : if(nspden==2)nitems=5 ! Total density, spin up, spin down, magnetization, zeta
1696 68932 : if(nspden==4)nitems=6 ! Total density, x, y, z, magnetization, zeta
1697 :
1698 137864 : ABI_MALLOC(value,(2,2,nitems))
1699 137864 : ABI_MALLOC(iindex,(2,2,nitems))
1700 206796 : ABI_MALLOC(array,(nfft))
1701 137864 : ABI_MALLOC(integrated,(nitems))
1702 :
1703 192119 : do iitems=1,nitems
1704 :
1705 : ! Copy the correct values into the array
1706 : ! First set of items : the density, for each spin component
1707 123187 : if(iitems<=nspden)then
1708 1780414225 : array(:)=rhor(:,iitems)
1709 : end if
1710 : ! Case nspden==2, some computation to be done
1711 123187 : if(nspden==2)then
1712 39525 : if(iitems==3)then ! Spin down
1713 207219197 : array(:)=rhor(:,1)-rhor(:,2)
1714 31620 : else if(iitems==4)then ! Magnetization
1715 207219197 : array(:)=2*rhor(:,2)-rhor(:,1)
1716 23715 : else if(iitems==5)then ! zeta = relative magnetization
1717 : ! Avoid 0/0: the limit of (x - y) / (x+ y) depends on the direction.
1718 207219197 : array(:)=zero
1719 207219197 : where (abs(rhor(:,1)) > tol12) array(:)=(2*rhor(:,2)-rhor(:,1))/rhor(:,1)
1720 : end if
1721 : ! Case nspden==4, some other computation to be done
1722 83662 : else if(nspden==4)then
1723 27162 : if(iitems==5)then ! Magnetization
1724 44154905 : array(:)=sqrt(rhor(:,2)**2+rhor(:,3)**2+rhor(:,4)**2)
1725 22635 : else if(iitems==6)then ! zeta = relative magnetization
1726 44154905 : array(:)=(sqrt(rhor(:,2)**2+rhor(:,3)**2+rhor(:,4)**2))/rhor(:,1)
1727 : end if
1728 : end if
1729 :
1730 : ! Zero all the absolute values that are lower than tol8, for portability reasons.
1731 2490381626 : do ifft = 1, nfft_
1732 2490381626 : if(abs(array(ifft))<tol8)array(ifft)=zero
1733 : end do
1734 :
1735 : ! DEBUG
1736 : ! write(std_out,*) ' iitems,array(1:2)=',iitems,array(1:2)
1737 : ! ENDDEBUG
1738 :
1739 369561 : do indsign=1,2 ! Find alternatively the maximum and the minimum
1740 246374 : iisign=3-2*indsign
1741 :
1742 246374 : if (nfft_ > 1) then
1743 : ! Initialize the two first values
1744 246374 : value1=array(istart + 1) ; value2=array(istart + 2)
1745 246374 : index1=1 ; index2=2
1746 :
1747 : ! Ordering, if needed
1748 246374 : if( iisign*(value2+tol12) > iisign*(value1)) then
1749 123187 : temp=value2 ; value2=value1 ; value1=temp
1750 123187 : index1=2 ; index2=1
1751 : end if
1752 :
1753 : ! Integration, if relevant
1754 246374 : if(present(ucvol).and. indsign==1)then
1755 123187 : integrated(iitems) = array(istart + 1)+array(istart + 2)
1756 : end if
1757 : else
1758 : value1 = zero; value2 = zero
1759 : index1 = 0; index2 = 0
1760 : end if
1761 :
1762 : ! DEBUG
1763 : ! write(std_out,*) ' value1,value2,index1,index2=',value1,value2,index1,index2
1764 : ! ENDDEBUG
1765 :
1766 : ! Loop over all points
1767 4980270504 : do ifft = 3, nfft_
1768 :
1769 4980024130 : temp=array(istart + ifft)
1770 4980024130 : if(present(ucvol).and. indsign==1)integrated(iitems) = integrated(iitems)+temp
1771 : ! Compares it to the second value
1772 4980270504 : if( iisign*(temp+tol12) > iisign*value2 ) then
1773 : ! Compare it to the first value
1774 210222812 : if( iisign*(temp+tol12) > iisign*value1 ) then
1775 207794951 : value2=value1 ; index2=index1
1776 207794951 : value1=temp ; index1=ifft
1777 : else
1778 : value2=temp ; index2=ifft
1779 : end if
1780 : end if
1781 :
1782 : end do ! ifft
1783 :
1784 246374 : value(1,indsign,iitems)=value1
1785 246374 : value(2,indsign,iitems)=value2
1786 246374 : iindex(1,indsign,iitems)=index1
1787 369561 : iindex(2,indsign,iitems)=index2
1788 :
1789 : ! DEBUG
1790 : ! write(std_out,*) ' it,v1,i1=',iitems, value1,index1
1791 : ! write(std_out,*) ' it,v2,i2=',iitems, value2,index2
1792 : ! ENDDEBUG
1793 :
1794 : end do ! indsign
1795 :
1796 192119 : if(present(ucvol))then
1797 123187 : nfftot=ngfft(1) * ngfft(2) * ngfft(3)
1798 123187 : integrated(iitems)=integrated(iitems)*ucvol/nfftot
1799 : end if
1800 :
1801 : ! Integrate the array
1802 : ! integrated(iitems)=zero
1803 : ! do ifft=1,nfft_
1804 : ! integrated(iitems) = integrated(iitems) + array(istart + ifft)
1805 : ! enddo
1806 : ! if(present(ucvol))integrated(iitems) = integrated(iitems)*ucvol/nfft_
1807 : ! write(std_err,*)present(ucvol)
1808 : ! if(present(ucvol))then
1809 : ! write(std_err,*)ucvol
1810 : ! endif
1811 :
1812 : end do ! iitems
1813 :
1814 68932 : ABI_FREE(array)
1815 :
1816 : !-------------------------------------------------------------------
1817 : !Enter section for FFT parallel case
1818 : !if(mpi_enreg%paral_kgb>1) spaceComm=mpi_enreg%comm_fft; reduce=.true.
1819 68932 : spaceComm=mpi_enreg%comm_fft; reduce=.false.
1820 68932 : if(mpi_enreg%nproc_fft>1) then
1821 4200 : spaceComm=mpi_enreg%comm_fft; reduce=.true.
1822 64732 : else if(mpi_enreg%nproc_wvl>1) then
1823 0 : spaceComm=mpi_enreg%comm_wvl; reduce=.true.
1824 : end if
1825 68932 : nproc=xmpi_comm_size(spaceComm)
1826 68932 : me=xmpi_comm_rank(spaceComm)
1827 :
1828 68932 : if (reduce) then
1829 :
1830 : ! Communicate all data to all processors with only two global communications
1831 16800 : ABI_MALLOC(value_fft,(5,nitems,nproc))
1832 16800 : ABI_MALLOC(index_fft,(2,2,nitems,nproc))
1833 72776 : value_fft(:,:,:)=zero
1834 82736 : index_fft(:,:,:,:)=0
1835 8972 : value_fft(1,:,me + 1)=value(1,1,:)
1836 8972 : value_fft(2,:,me + 1)=value(2,1,:)
1837 8972 : value_fft(3,:,me + 1)=value(1,2,:)
1838 8972 : value_fft(4,:,me + 1)=value(2,2,:)
1839 8972 : if(present(ucvol))value_fft(5,:,me + 1)=integrated(:)
1840 37604 : index_fft(:,:,:,me + 1)=iindex(:,:,:)
1841 4200 : call xmpi_sum(value_fft,spaceComm,ierr)
1842 4200 : call xmpi_sum(index_fft,spaceComm,ierr)
1843 :
1844 : ! Determine the global optimum and second optimum for each item
1845 : ! Also, the integrated quantities, if relevant.
1846 8972 : do iitems=1,nitems
1847 :
1848 14732 : if(present(ucvol))integrated(iitems)=sum(value_fft(5,iitems,1:nproc))
1849 :
1850 18516 : do indsign=1,2 ! Find alternatively the maximum and the minimum
1851 9544 : iisign=3-2*indsign
1852 :
1853 : ! Initialisation
1854 9544 : value1=value_fft(2*indsign-1,iitems,1)
1855 9544 : value2=value_fft(2*indsign ,iitems,1)
1856 9544 : index1=index_fft(1,indsign,iitems,1)
1857 9544 : index2=index_fft(2,indsign,iitems,1)
1858 :
1859 : ! Loop
1860 29464 : do iproc=1, nproc, 1
1861 69304 : do ii=1,2
1862 59760 : if(iproc>1 .or. ii==2)then
1863 :
1864 30296 : temp=value_fft(ii+2*(indsign-1),iitems,iproc)
1865 : ! Compares it to the second value
1866 30296 : if( iisign*(temp+tol12) > iisign*value2 ) then
1867 : ! Compare it to the first value
1868 13040 : if( iisign*(temp+tol12) > iisign*value1 ) then
1869 6384 : value2=value1 ; index2=index1
1870 6384 : value1=temp ; index1=index_fft(ii,indsign,iitems,iproc)
1871 : else
1872 6656 : value2=temp ; index2=index_fft(ii,indsign,iitems,iproc)
1873 : end if
1874 : end if
1875 :
1876 : end if ! if(iproc>1 .or. ii==2)
1877 : end do ! ii
1878 : end do ! iproc
1879 :
1880 9544 : value(1,indsign,iitems)=value1
1881 9544 : value(2,indsign,iitems)=value2
1882 9544 : iindex(1,indsign,iitems)=index1
1883 14316 : iindex(2,indsign,iitems)=index2
1884 :
1885 : end do ! iisign
1886 : end do ! iitems
1887 :
1888 4200 : ABI_FREE(value_fft)
1889 8400 : ABI_FREE(index_fft)
1890 :
1891 : end if !if(reduce)
1892 :
1893 : !-------------------------------------------------------------------
1894 :
1895 : !Determines the reduced coordinates of the min and max for each item
1896 137864 : ABI_MALLOC(coord,(3,2,2,nitems))
1897 192119 : do iitems=1,nitems
1898 438493 : do indsign=1,2
1899 862309 : do ii=1,2
1900 492748 : index1=iindex(ii,indsign,iitems)
1901 492748 : i3=(index1-1)/n1/n2
1902 492748 : i2=(index1-1-i3*n1*n2)/n1
1903 492748 : i1=index1-1-i3*n1*n2-i2*n1
1904 492748 : coord(1,ii,indsign,iitems)=dble(i1)/dble(n1)+tol12
1905 492748 : coord(2,ii,indsign,iitems)=dble(i2)/dble(n2)+tol12
1906 739122 : coord(3,ii,indsign,iitems)=dble(i3)/dble(n3)+tol12
1907 : ! DEBUG
1908 : ! write(std_out,*)' ii,indsign,iitems,coord(1:3)=',ii,indsign,iitems,coord(:,ii,indsign,iitems)
1909 : ! write(std_out,*)' value ', value(ii, indsign, iitems)
1910 : ! ENDDEBUG
1911 : end do
1912 : end do
1913 : end do
1914 :
1915 : !-------------------------------------------------------------------------
1916 : !Output
1917 68932 : if (mpi_enreg%paral_kgb==0.or.mpi_enreg%me_fft==0) then
1918 : if(.true.)then
1919 187523 : do iitems=1,nitems
1920 :
1921 120745 : if(ioptrhor==4 .and. iitems>2)exit
1922 :
1923 119133 : select case (ioptrhor)
1924 : case(0)
1925 :
1926 119133 : if(iitems==1) write(message,'(a)')' Total charge density [el/Bohr^3]'
1927 119133 : if(nspden==2)then
1928 38790 : if(iitems==2) write(message,'(a)')' Spin up density [el/Bohr^3]'
1929 38790 : if(iitems==3) write(message,'(a)')' Spin down density [el/Bohr^3]'
1930 38790 : if(iitems==4) write(message,'(a)')' Magnetization (spin up - spin down) [el/Bohr^3]'
1931 38790 : if(iitems==5) write(message,'(a)')' Relative magnetization (=zeta, between -1 and 1) '
1932 80343 : else if(nspden==4)then
1933 26646 : if(iitems==2) write(message,'(a)')' x component of magnetization [el/Bohr^3]'
1934 26646 : if(iitems==3) write(message,'(a)')' y component of magnetization [el/Bohr^3]'
1935 26646 : if(iitems==4) write(message,'(a)')' z component of magnetization [el/Bohr^3]'
1936 26646 : if(iitems==5) write(message,'(a)')' Magnetization (absolute value) [el/Bohr^3]'
1937 26646 : if(iitems==6) write(message,'(a)')' Relative magnetization (=zeta, between -1 and 1) '
1938 : end if
1939 :
1940 : case(1)
1941 :
1942 1577 : if(iitems==1) write(message,'(a)')' Total kinetic energy density [Ha/Bohr^3]'
1943 1577 : if(nspden==2)then
1944 535 : if(iitems==2) write(message,'(a)')' Spin up density [Ha/Bohr^3]'
1945 535 : if(iitems==3) write(message,'(a)')' Spin down density [Ha/Bohr^3]'
1946 535 : if(iitems==4) write(message,'(a)')' Magnetization (spin up - spin down) [Ha/Bohr^3]'
1947 535 : if(iitems==5) write(message,'(a)')' Relative magnetization (=zeta, between -1 and 1) '
1948 1042 : else if(nspden==4)then
1949 336 : if(iitems==2) write(message,'(a)')' x component of magnetization [Ha/Bohr^3]'
1950 336 : if(iitems==3) write(message,'(a)')' y component of magnetization [Ha/Bohr^3]'
1951 336 : if(iitems==4) write(message,'(a)')' z component of magnetization [Ha/Bohr^3]'
1952 336 : if(iitems==5) write(message,'(a)')' Magnetization (absolute value) [Ha/Bohr^3]'
1953 336 : if(iitems==6) write(message,'(a)')' Relative magnetization (=zeta, between -1 and 1) '
1954 : end if
1955 :
1956 : case(2)
1957 :
1958 21 : if(iitems==1) write(message,'(a)')' Gradient of the electronic density [el/Bohr^4]'
1959 21 : if(nspden==2)then
1960 15 : if(iitems==2) write(message,'(a)')' Spin up density [el/Bohr^4]'
1961 15 : if(iitems==3) write(message,'(a)')' Spin down density [el/Bohr^4]'
1962 15 : if(iitems==4) write(message,'(a)')' Magnetization (spin up - spin down) [el/Bohr^4]'
1963 15 : if(iitems==5) write(message,'(a)')' Relative magnetization (=zeta, between -1 and 1) '
1964 6 : else if(nspden==4)then
1965 0 : if(iitems==2) write(message,'(a)')' x component of magnetization [el/Bohr^4]'
1966 0 : if(iitems==3) write(message,'(a)')' y component of magnetization [el/Bohr^4]'
1967 0 : if(iitems==4) write(message,'(a)')' z component of magnetization [el/Bohr^4]'
1968 0 : if(iitems==5) write(message,'(a)')' Magnetization (absolute value) [el/Bohr^4]'
1969 0 : if(iitems==6) write(message,'(a)')' Relative magnetization (=zeta, between -1 and 1) '
1970 : end if
1971 :
1972 : case(3)
1973 :
1974 6 : if(iitems==1) write(message,'(a)')' Laplacian of the electronic density [el/Bohr^5]'
1975 6 : if(nspden==2)then
1976 5 : if(iitems==2) write(message,'(a)')' Spin up density [el/Bohr^5]'
1977 5 : if(iitems==3) write(message,'(a)')' Spin down density [el/Bohr^5]'
1978 5 : if(iitems==4) write(message,'(a)')' Magnetization (spin up - spin down) [el/Bohr^5]'
1979 5 : if(iitems==5) write(message,'(a)')' Relative magnetization (=zeta, between -1 and 1) '
1980 1 : else if(nspden==4)then
1981 0 : if(iitems==2) write(message,'(a)')' x component of magnetization [el/Bohr^5]'
1982 0 : if(iitems==3) write(message,'(a)')' y component of magnetization [el/Bohr^5]'
1983 0 : if(iitems==4) write(message,'(a)')' z component of magnetization [el/Bohr^5]'
1984 0 : if(iitems==5) write(message,'(a)')' Magnetization (absolute value) [el/Bohr^5]'
1985 0 : if(iitems==6) write(message,'(a)')' Relative magnetization (=zeta, between -1 and 1) '
1986 : end if
1987 :
1988 : case(4)
1989 :
1990 6 : if(iitems==1) write(message,'(a)')' Electron Localization Function (ELF) [min:0;max:1]'
1991 120749 : if(nspden==2)then
1992 4 : if(iitems==2) write(message,'(a)')' Spin up ELF [min:0;max:1]'
1993 : ! if(iitems==3) write(message,'(a)')' Spin down ELF [min:0;max:1]'
1994 : ! if(iitems==4) write(message,'(a)')' Magnetization (spin up - spin down) [el/Bohr^4]'
1995 : ! if(iitems==5) write(message,'(a)')' Relative magnetization (=zeta, between -1 and 1) '
1996 : else if(nspden==4)then
1997 : ! if(iitems==2) write(message,'(a)')' x component of magnetization [el/Bohr^4]'
1998 : ! if(iitems==3) write(message,'(a)')' y component of magnetization [el/Bohr^4]'
1999 : ! if(iitems==4) write(message,'(a)')' z component of magnetization [el/Bohr^4]'
2000 : ! if(iitems==5) write(message,'(a)')' Magnetization (spin up - spin down) [el/Bohr^4]'
2001 : ! if(iitems==6) write(message,'(a)')' Relative magnetization (=zeta, between -1 and 1) '
2002 : end if
2003 : end select
2004 :
2005 120743 : call wrtout(iout,message,'COLL')
2006 :
2007 120743 : write(message,'(a,es13.4,a,3f10.4)') ') Maximum= ',&
2008 241486 : & value(1,1,iitems),' at reduced coord.',coord(:,1,1,iitems)
2009 120743 : call wrtout(iout,message,'COLL')
2010 120743 : if(option==2)then
2011 9367 : write(message,'(a,es13.4,a,3f10.4)') ')Next maximum= ',&
2012 18734 : & value(2,1,iitems),' at reduced coord.',coord(:,2,1,iitems)
2013 9367 : call wrtout(iout,message,'COLL')
2014 : end if
2015 120743 : write(message,'(a,es13.4,a,3f10.4)') ') Minimum= ',&
2016 241486 : & value(1,2,iitems),' at reduced coord.',coord(:,1,2,iitems)
2017 120743 : call wrtout(iout,message,'COLL')
2018 120743 : if(option==2)then
2019 9367 : write(message,'(a,es13.4,a,3f10.4)') ')Next minimum= ',&
2020 18734 : & value(2,2,iitems),' at reduced coord.',coord(:,2,2,iitems)
2021 9367 : call wrtout(iout,message,'COLL')
2022 : end if
2023 187523 : if(present(ucvol))then
2024 120743 : if(.not.(nspden==2.and.iitems==5) .and. .not.(nspden==4.and.iitems==6))then
2025 108377 : if(abs(integrated(iitems))<tol10)integrated(iitems)=zero
2026 108377 : write(message,'(a,es13.4)')' Integrated= ',integrated(iitems)
2027 108377 : call wrtout(iout,message,'COLL')
2028 : end if
2029 : end if
2030 :
2031 : end do ! iitems
2032 : end if
2033 :
2034 : if(.false.)then
2035 :
2036 : select case(optrhor)
2037 : case(0)
2038 : write(txt1_in_mssg, '(a)')" Min el dens="
2039 : write(txt2_in_mssg, '(a)')" el/bohr^3 at reduced coord."
2040 : write(txt3_in_mssg, '(a)')" Max el dens="
2041 : case(1)
2042 : write(txt1_in_mssg, '(a)')" Min kin energy dens="
2043 : write(txt2_in_mssg, '(a)')" bohr^(-5) at reduced coord."
2044 : write(txt3_in_mssg, '(a)')" Max kin energy dens="
2045 : end select
2046 :
2047 : write(message, '(a,a,1p,e12.4,a,0p,3f8.4)' ) ch10,&
2048 : & trim(txt1_in_mssg),value(1,2,1),&
2049 : & trim(txt2_in_mssg),coord(:,1,2,1)
2050 : call wrtout(iout,message,'COLL')
2051 : if(option==2)then
2052 : write(message, '(a,1p,e12.4,a,0p,3f8.4)' ) &
2053 : & ', next min=',value(2,2,1),&
2054 : & trim(txt2_in_mssg),coord(:,2,2,1)
2055 : call wrtout(iout,message,'COLL')
2056 : end if
2057 : write(message, '(a,1p,e12.4,a,0p,3f8.4)' )&
2058 : & trim(txt3_in_mssg),value(1,1,1),&
2059 : & trim(txt2_in_mssg),coord(:,1,1,1)
2060 : call wrtout(iout,message,'COLL')
2061 : if(option==2)then
2062 : write(message, '(a,1p,e12.4,a,0p,3f8.4)' )&
2063 : & ', next max=',value(2,1,1),&
2064 : & trim(txt2_in_mssg),coord(:,2,1,1)
2065 : call wrtout(iout,message,'COLL')
2066 : end if
2067 :
2068 : if(nspden>=2)then
2069 : write(message, '(a,a,1p,e12.4,a,0p,3f8.4)' ) ch10,&
2070 : & ',Min spin pol zeta=',value(1,2,4+nspden/2),&
2071 : & ' at reduced coord.',coord(:,1,2,4+nspden/2)
2072 : call wrtout(iout,message,'COLL')
2073 : if(option==2)then
2074 : write(message, '(a,1p,e12.4,a,0p,3f8.4)' )&
2075 : & ', next min=',value(2,2,4+nspden/2),&
2076 : & ' at reduced coord.',coord(:,2,2,4+nspden/2)
2077 : call wrtout(iout,message,'COLL')
2078 : end if
2079 : write(message, '(a,1p,e12.4,a,0p,3f8.4)' )&
2080 : & ',Max spin pol zeta=',value(1,1,4+nspden/2),&
2081 : & ' at reduced coord.',coord(:,1,1,4+nspden/2)
2082 : call wrtout(iout,message,'COLL')
2083 : if(option==2)then
2084 : write(message, '(a,1p,e12.4,a,0p,3f8.4)' )&
2085 : & ', next max=',value(2,1,4+nspden/2),&
2086 : & ' at reduced coord.',coord(:,2,1,4+nspden/2)
2087 : call wrtout(iout,message,'COLL')
2088 : end if
2089 : end if ! nspden
2090 :
2091 : end if ! second section always true
2092 :
2093 : if(nspden==2 .and. .false.)then
2094 : write(message,'(a)')&
2095 : & ' Position in reduced coord. ( x y z )'
2096 : call wrtout(iout,message,'COLL')
2097 : write(message,'(a,es13.4,a,3f10.4)')' Minimum (Total el-den) : [el/Bohr^3]',&
2098 : & rhomn1(1),' at',ri_rhomn1(1,1),ri_rhomn1(2,1),ri_rhomn1(3,1)
2099 : call wrtout(iout,message,'COLL')
2100 : write(message,'(a,es13.4,a,3f10.4)')' Minimum (Spin-up den) : [el/Bohr^3]',&
2101 : & rhomn1(2),' at',ri_rhomn1(1,2),ri_rhomn1(2,2),ri_rhomn1(3,2)
2102 : call wrtout(iout,message,'COLL')
2103 : write(message,'(a,es13.4,a,3f10.4)')' Minimum (Spin-down den) : [el/Bohr^3]',&
2104 : & zetmn1(1),' at',ri_zetmn1(1,1),ri_zetmn1(2,1),ri_zetmn1(3,1)
2105 : call wrtout(iout,message,'COLL')
2106 : write(message,'(a,es13.4,a,3f10.4)')' Minimum (Spin pol zeta) : [m/|m|] ',&
2107 : & zetmn1(2),' at',ri_zetmn1(1,2),ri_zetmn1(2,2),ri_zetmn1(3,2)
2108 : call wrtout(iout,message,'COLL')
2109 : if(option==2)then
2110 : write(message,'(a,es13.4,a,3f10.4)')' Next minimum (Total el-den) : [el/Bohr^3]',&
2111 : & rhomn2(1),' at',ri_rhomn2(1,1),ri_rhomn2(2,1),ri_rhomn2(3,1)
2112 : call wrtout(iout,message,'COLL')
2113 : write(message,'(a,es13.4,a,3f10.4)')' Next minimum (Spin-up den) : [el/Bohr^3]',&
2114 : & rhomn2(2),' at',ri_rhomn2(1,2),ri_rhomn2(2,2),ri_rhomn2(3,2)
2115 : call wrtout(iout,message,'COLL')
2116 : write(message,'(a,es13.4,a,3f10.4)')' Next minimum (Spin-down den) : [el/Bohr^3]',&
2117 : & zetmn2(1),' at',ri_zetmn2(1,1),ri_zetmn2(2,1),ri_zetmn2(3,1)
2118 : call wrtout(iout,message,'COLL')
2119 : write(message,'(a,es13.4,a,3f10.4)')' Next minimum (Spin pol zeta) : [m/|m|] ',&
2120 : & zetmn2(2),' at',ri_zetmn2(1,2),ri_zetmn2(2,2),ri_zetmn2(3,2)
2121 : call wrtout(iout,message,'COLL')
2122 : end if
2123 : write(message,*)' '
2124 : call wrtout(iout,message,'COLL')
2125 : write(message,'(a,es13.4,a,3f10.4)')' Maximum (Total el-den) : [el/Bohr^3]',&
2126 : & rhomx1(1),' at',ri_rhomx1(1,1),ri_rhomx1(2,1),ri_rhomx1(3,1)
2127 : call wrtout(iout,message,'COLL')
2128 : write(message,'(a,es13.4,a,3f10.4)')' Maximum (Spin-up den) : [el/Bohr^3]',&
2129 : & rhomx1(2),' at',ri_rhomx1(1,2),ri_rhomx1(2,2),ri_rhomx1(3,2)
2130 : call wrtout(iout,message,'COLL')
2131 : write(message,'(a,es13.4,a,3f10.4)')' Maximum (Spin-down den) : [el/Bohr^3]',&
2132 : & zetmx1(1),' at',ri_zetmx1(1,1),ri_zetmx1(2,1),ri_zetmx1(3,1)
2133 : call wrtout(iout,message,'COLL')
2134 : write(message,'(a,es13.4,a,3f10.4)')' Maximum (Spin pol zeta) : [m/|m|] ',&
2135 : & zetmx1(2),' at',ri_zetmx1(1,2),ri_zetmx1(2,2),ri_zetmx1(3,2)
2136 : call wrtout(iout,message,'COLL')
2137 : if(option==2)then
2138 : write(message,'(a,es13.4,a,3f10.4)')' Next maximum (Total el-den) : [el/Bohr^3]',&
2139 : & rhomx2(1),' at',ri_rhomx2(1,1),ri_rhomx2(2,1),ri_rhomx2(3,1)
2140 : call wrtout(iout,message,'COLL')
2141 : write(message,'(a,es13.4,a,3f10.4)')' Next maximum (Spin-up den) : [el/Bohr^3]',&
2142 : & rhomx2(2),' at',ri_rhomx2(1,2),ri_rhomx2(2,2),ri_rhomx2(3,2)
2143 : call wrtout(iout,message,'COLL')
2144 : write(message,'(a,es13.4,a,3f10.4)')' Next maximum (Spin-down den) : [el/Bohr^3]',&
2145 : & zetmx2(1),' at',ri_zetmx2(1,1),ri_zetmx2(2,1),ri_zetmx2(3,1)
2146 : call wrtout(iout,message,'COLL')
2147 : write(message,'(a,es13.4,a,3f10.4)')' Next maximum (Spin pol zeta) : [m/|m|] ',&
2148 : & zetmx2(2),' at',ri_zetmx2(1,2),ri_zetmx2(2,2),ri_zetmx2(3,2)
2149 : call wrtout(iout,message,'COLL')
2150 : end if
2151 : end if
2152 :
2153 : if(nspden==4 .and. .false.)then
2154 : write(message,'(a)')&
2155 : & ' Position in reduced coord. ( x y z )'
2156 : call wrtout(iout,message,'COLL')
2157 : write(message,'(a,es13.4,a,3f10.4)')' Minimum (Total el-den) : [el/Bohr^3]',&
2158 : & rhomn1(1),' at',ri_rhomn1(1,1),ri_rhomn1(2,1),ri_rhomn1(3,1)
2159 : call wrtout(iout,message,'COLL')
2160 : write(message,'(a,es13.4,a,3f10.4)')' Minimum (Magnetizat.-x) : [m/|m|] ',&
2161 : & rhomn1(2),' at',ri_rhomn1(1,2),ri_rhomn1(2,2),ri_rhomn1(3,2)
2162 : call wrtout(iout,message,'COLL')
2163 : write(message,'(a,es13.4,a,3f10.4)')' Minimum (Magnetizat.-y) : [m/|m|] ',&
2164 : & rhomn1(3),' at',ri_rhomn1(1,3),ri_rhomn1(2,3),ri_rhomn1(3,3)
2165 : call wrtout(iout,message,'COLL')
2166 : write(message,'(a,es13.4,a,3f10.4)')' Minimum (Magnetizat.-z) : [m/|m|] ',&
2167 : & rhomn1(4),' at',ri_rhomn1(1,4),ri_rhomn1(2,4),ri_rhomn1(3,4)
2168 : call wrtout(iout,message,'COLL')
2169 : write(message,'(a,es13.4,a,3f10.4)')' Minimum (Spin pol zeta) : [m/|m|] ',&
2170 : & zetmn1(1),' at',ri_zetmn1(1,1),ri_zetmn1(2,1),ri_zetmn1(3,1)
2171 : call wrtout(iout,message,'COLL')
2172 : if(option==2)then
2173 : write(message,'(a,es13.4,a,3f10.4)')' Next-Minimum (Total el-den) : [el/Bohr^3]',&
2174 : & rhomn2(1),' at',ri_rhomn2(1,1),ri_rhomn2(2,1),ri_rhomn2(3,1)
2175 : call wrtout(iout,message,'COLL')
2176 : write(message,'(a,es13.4,a,3f10.4)')' Next-Minimum (Magnetizat.-x) : [m/|m|] ',&
2177 : & rhomn2(2),' at',ri_rhomn2(1,2),ri_rhomn2(2,2),ri_rhomn2(3,2)
2178 : call wrtout(iout,message,'COLL')
2179 : write(message,'(a,es13.4,a,3f10.4)')' Next-Minimum (Magnetizat.-y) : [m/|m|] ',&
2180 : & rhomn2(3),' at',ri_rhomn2(1,3),ri_rhomn2(2,3),ri_rhomn2(3,3)
2181 : call wrtout(iout,message,'COLL')
2182 : write(message,'(a,es13.4,a,3f10.4)')' Next-Minimum (Magnetizat.-z) : [m/|m|] ',&
2183 : & rhomn2(4),' at',ri_rhomn2(1,4),ri_rhomn2(2,4),ri_rhomn2(3,4)
2184 : call wrtout(iout,message,'COLL')
2185 : write(message,'(a,es13.4,a,3f10.4)')' Next-Minimum (Spin pol zeta) : [m/|m|] ',&
2186 : & zetmn2(1),' at',ri_zetmn2(1,1),ri_zetmn2(2,1),ri_zetmn2(3,1)
2187 : call wrtout(iout,message,'COLL')
2188 : end if
2189 : write(message,*)' '
2190 : call wrtout(iout,message,'COLL')
2191 : write(message,'(a,es13.4,a,3f10.4)')' Maximum (Total el-den) : [el/Bohr^3]',&
2192 : & rhomx1(1),' at',ri_rhomx1(1,1),ri_rhomx1(2,1),ri_rhomx1(3,1)
2193 : call wrtout(iout,message,'COLL')
2194 : write(message,'(a,es13.4,a,3f10.4)')' Maximum (Magnetizat.-x) : [m/|m|] ',&
2195 : & rhomx1(2),' at',ri_rhomx1(1,2),ri_rhomx1(2,2),ri_rhomx1(3,2)
2196 : call wrtout(iout,message,'COLL')
2197 : write(message,'(a,es13.4,a,3f10.4)')' Maximum (Magnetizat.-y) : [m/|m|] ',&
2198 : & rhomx1(3),' at',ri_rhomx1(1,3),ri_rhomx1(2,3),ri_rhomx1(3,3)
2199 : call wrtout(iout,message,'COLL')
2200 : write(message,'(a,es13.4,a,3f10.4)')' Maximum (Magnetizat.-z) : [m/|m|] ',&
2201 : & rhomx1(4),' at',ri_rhomx1(1,4),ri_rhomx1(2,4),ri_rhomx1(3,4)
2202 : call wrtout(iout,message,'COLL')
2203 : write(message,'(a,es13.4,a,3f10.4)')' Maximum (Spin pol zeta) : [m/|m|] ',&
2204 : & zetmx1(1),' at',ri_zetmx1(1,1),ri_zetmx1(2,1),ri_zetmx1(3,1)
2205 : call wrtout(iout,message,'COLL')
2206 : if(option==2)then
2207 : write(message,'(a,es13.4,a,3f10.4)')' Next-Maximum (Total el-den) : [el/Bohr^3]',&
2208 : & rhomx2(1),' at',ri_rhomx2(1,1),ri_rhomx2(2,1),ri_rhomx2(3,1)
2209 : call wrtout(iout,message,'COLL')
2210 : write(message,'(a,es13.4,a,3f10.4)')' Next-Maximum (Magnetizat.-x) : [m/|m|] ',&
2211 : & rhomx2(2),' at',ri_rhomx2(1,2),ri_rhomx2(2,2),ri_rhomx2(3,2)
2212 : call wrtout(iout,message,'COLL')
2213 : write(message,'(a,es13.4,a,3f10.4)')' Next-Maximum (Magnetizat.-y) : [m/|m|] ',&
2214 : & rhomx2(3),' at',ri_rhomx2(1,3),ri_rhomx2(2,3),ri_rhomx2(3,3)
2215 : call wrtout(iout,message,'COLL')
2216 : write(message,'(a,es13.4,a,3f10.4)')' Next-Maximum (Magnetizat.-z) : [m/|m|] ',&
2217 : & rhomx2(4),' at',ri_rhomx2(1,4),ri_rhomx2(2,4),ri_rhomx2(3,4)
2218 : call wrtout(iout,message,'COLL')
2219 : write(message,'(a,es13.4,a,3f10.4)')' Next-Maximum (Spin pol zeta) : [m/|m|] ',&
2220 : & zetmx2(1),' at',ri_zetmx2(1,1),ri_zetmx2(2,1),ri_zetmx2(3,1)
2221 : call wrtout(iout,message,'COLL')
2222 : end if
2223 : end if
2224 : end if
2225 :
2226 68932 : ABI_FREE(coord)
2227 68932 : ABI_FREE(value)
2228 68932 : ABI_FREE(iindex)
2229 68932 : ABI_FREE(integrated)
2230 :
2231 68932 : end subroutine prtrhomxmn
2232 : !!***
2233 :
2234 : !!****f* m_mkrho/read_atomden
2235 : !! NAME
2236 : !! read_atomden
2237 : !!
2238 : !! FUNCTION
2239 : !!
2240 : !! INPUTS
2241 : !! natom : number of atoms in cell
2242 : !! nfft=(effective) number of FFT grid points (for this processor) - fine grid
2243 : !! ngfft(18)=contain all needed information about 3D FFT,
2244 : !! nspden : number of spin densities
2245 : !! ntypat : number of types of atoms in the cell
2246 : !! typat(natom) : list of atom types
2247 : !!
2248 : !! OUTPUT
2249 : !! rhor_atm(nfft,nspden) : full electron density on the (fine) grid
2250 : !!
2251 : !! SOURCE
2252 :
2253 0 : subroutine read_atomden(MPI_enreg,natom,nfft,ngfft,nspden,ntypat, &
2254 0 : & rhor_atm,typat,rprimd,xred,prtvol,file_prefix)
2255 :
2256 : !Arguments ------------------------------------
2257 : !scalars
2258 : integer,intent(in) :: natom,nfft,nspden,ntypat,prtvol
2259 : !arrays
2260 : type(MPI_type),intent(in) :: MPI_enreg
2261 : integer,intent(in) :: ngfft(18),typat(natom)
2262 : real(dp), intent(in) :: rprimd(3,3),xred(3,natom)
2263 : real(dp),intent(inout) :: rhor_atm(nfft,nspden)
2264 : character(len=7), intent(in) :: file_prefix
2265 :
2266 : !Local variables-------------------------------
2267 : !scalars
2268 : character(len=500) :: message
2269 : character(len=120) :: filename
2270 : character(len=7) :: calctype='replace'
2271 : integer :: igrid,i,i1,i2,i3,io_err,itypat,unt
2272 : integer :: natomgrmax,nlines,ngrid,n1,n2,n3
2273 : real(dp) :: difx,dify,difz,ucvol!,norm
2274 : !arrays
2275 0 : integer :: natomgr(ntypat)
2276 : real(dp) :: a(3),b(3),c(3)
2277 : real(dp) :: gmet(3,3),gprimd(3,3),rmet(3,3)
2278 0 : real(dp),allocatable :: atomrgrid(:,:),r_vec_grid(:,:),density(:,:)
2279 0 : real(dp),allocatable :: rho(:)
2280 :
2281 : ! ************************************************************************
2282 :
2283 : !Initialise various variables
2284 0 : ngrid = nfft
2285 0 : a(:) = rprimd(:,1)
2286 0 : b(:) = rprimd(:,2)
2287 0 : c(:) = rprimd(:,3)
2288 0 : ABI_MALLOC(rho,(ngrid))
2289 0 : if (nspden/=1) then
2290 0 : ABI_ERROR('read_atomden: Only nspden=1 allowed.')
2291 : end if
2292 0 : rho = rhor_atm(:,1)
2293 0 : gmet=zero;gprimd=zero;rmet=zero;ucvol=zero
2294 :
2295 :
2296 : !Calculate the r vector (reduced coord.) of the fine gridpoints
2297 0 : ABI_MALLOC(r_vec_grid,(3,ngrid))
2298 0 : igrid = 0
2299 0 : n1 = ngfft(1)
2300 0 : n2 = ngfft(2)
2301 0 : n3 = ngfft(3)
2302 0 : do i3=0,n3-1
2303 0 : difz=dble(i3)/dble(n3)
2304 0 : do i2=0,n2-1
2305 0 : dify=dble(i2)/dble(n2)
2306 0 : do i1=0,n1-1
2307 0 : difx=dble(i1)/dble(n1)
2308 0 : igrid = igrid + 1
2309 0 : r_vec_grid(1,igrid)=difx*rprimd(1,1)+dify*rprimd(1,2)+difz*rprimd(1,3)
2310 0 : r_vec_grid(2,igrid)=difx*rprimd(2,1)+dify*rprimd(2,2)+difz*rprimd(2,3)
2311 0 : r_vec_grid(3,igrid)=difx*rprimd(3,1)+dify*rprimd(3,2)+difz*rprimd(3,3)
2312 : end do
2313 : end do
2314 : end do
2315 0 : if (igrid/=ngrid) then
2316 0 : ABI_ERROR('read_atomden: igrid not equal to ngrid')
2317 : end if
2318 :
2319 : !Read in atomic density data for each atom type
2320 : !first check how many datapoints are in each file
2321 0 : do itypat=1,ntypat
2322 0 : filename='';io_err=0;
2323 0 : if (itypat>0) write(filename,'(a,a,i1,a)') trim(file_prefix), '_density_atom_type',itypat,'.dat'
2324 0 : if (itypat>10) write(filename,'(a,a,i2,a)') trim(file_prefix), '_density_atom_type',itypat,'.dat'
2325 0 : if (open_file(filename, message, newunit=unt, status='old',action='read') /= 0) then
2326 0 : write(std_out,*) 'ERROR in read_atomden: Could not open file: ',filename
2327 0 : write(std_out,*) ' Current implementation requires this file to be present'
2328 0 : write(std_out,*) ' for each type of atom.'
2329 0 : write(std_out,*)trim(message)
2330 0 : ABI_ERROR("Cannot continue")
2331 : end if
2332 : ! Check number of lines in file
2333 0 : nlines = 1;io_err=0;
2334 0 : do
2335 0 : read(unt,*,iostat=io_err)
2336 0 : if (io_err<0) exit
2337 0 : nlines = nlines + 1
2338 : end do
2339 0 : close(unt)
2340 0 : natomgr(itypat) = nlines - 2
2341 : end do ! Atom type
2342 : !Allocate arrays and read in data
2343 0 : natomgrmax = maxval(natomgr)
2344 0 : ABI_MALLOC(atomrgrid,(natomgrmax,ntypat))
2345 0 : ABI_MALLOC(density,(natomgrmax,ntypat))
2346 0 : atomrgrid = zero ; density = zero
2347 0 : do itypat=1,ntypat
2348 0 : filename='';io_err=0;
2349 0 : if (itypat>0) write(filename,'(a,a,i1,a)') trim(file_prefix), '_density_atom_type',itypat,'.dat'
2350 0 : if (itypat>10) write(filename,'(a,a,i2,a)') trim(file_prefix), '_density_atom_type',itypat,'.dat'
2351 0 : if (open_file(filename,message,newunit=unt,status='old',action='read') /= 0) then
2352 0 : ABI_ERROR(message)
2353 : end if
2354 0 : read(unt,*) ! Skip comment line
2355 0 : do i=1,natomgr(itypat)
2356 0 : read(unt,*) atomrgrid(i,itypat),density(i,itypat)
2357 : end do
2358 0 : close(unt)
2359 0 : if (atomrgrid(1,itypat)/=zero) then
2360 0 : write(std_out,*) 'ERROR in read_atomden, in file: ',filename
2361 0 : write(std_out,*) ' First gridpoint has to be the origin.'
2362 0 : ABI_ERROR("Cannot continue")
2363 : end if
2364 : end do ! Atom type
2365 :
2366 : !write(std_out,*) '*** --- In read_atomden before call--- ***'
2367 : !write(std_out,*) ' calctype:',calctype,' natom:',natom
2368 : !write(std_out,*) ' ntypat:',ntypat,' typat:',typat
2369 : !write(std_out,*) ' ngrid:',ngrid
2370 : !write(std_out,*) ' a:',a
2371 : !write(std_out,*) ' b:',b
2372 : !write(std_out,*) ' c:',c
2373 : !write(std_out,*) ' xred:',xred
2374 : !write(std_out,*) ' natomgr:',natomgr
2375 : !write(std_out,*) 'natomgrmax:',natomgrmax
2376 : !write(std_out,*) ' atomrgrid:',atomrgrid
2377 : !write(std_out,*) ' density:',density
2378 : !write(std_out,*) 'r_vec_grid:'
2379 : !write(std_out,*) r_vec_grid
2380 :
2381 : !Call atomden
2382 : call atomden(MPI_enreg,natom,ntypat,typat,ngrid,r_vec_grid,rho,a,b,c,xred, &
2383 0 : & natomgr,natomgrmax,atomrgrid,density,prtvol,calctype)
2384 :
2385 : !if (prtvol>9) then ! calculate norm
2386 : !call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
2387 : !norm = SUM(rho(:))*ucvol/dble(n1*n2*n3)
2388 : !write(message,'(a,F8.4)') ' In read_atomden - NORM OF DENSITY: ',norm
2389 : !call wrtout(std_out,message,'COLL')
2390 : !end if
2391 :
2392 0 : rhor_atm(:,1) = rho
2393 :
2394 0 : ABI_SFREE(atomrgrid)
2395 0 : ABI_SFREE(density)
2396 0 : ABI_SFREE(r_vec_grid)
2397 0 : ABI_SFREE(rho)
2398 :
2399 0 : end subroutine read_atomden
2400 : !!***
2401 :
2402 : !!****f* m_mkrho/atomden
2403 : !! NAME
2404 : !! atomden
2405 : !!
2406 : !! FUNCTION
2407 : !! Construct atomic proto-bulk density (i.e. the superposed density
2408 : !! from neutral, isolated atoms at the bulk atomic positions).
2409 : !! This is useful if one wants to construct the bonding density:
2410 : !!
2411 : !! rho^{bnd} = rho^{bulk}(r)
2412 : !! - \sum_{\alpha}\rho^{atm}_{\alpha}(r-R_{\alpha})
2413 : !!
2414 : !! Where rho^{bulk} is the bulk density, rho^{atm} the atomic density
2415 : !! and the index \alpha sums over all atoms. the R_{\alpha} are the
2416 : !! atomic positions in the bulk. This routine calculates the sum over
2417 : !! rho^{atm}_{\alpha}(r-R_{\alpha}) on a grid.
2418 : !!
2419 : !! Units are atomic.
2420 : !!
2421 : !! INPUTS
2422 : !! calctype : type of calculation
2423 : !! 'replace' zero the input/output density array
2424 : !! 'add' add to the input/output density array
2425 : !! natom : number of atoms in cell
2426 : !! ntypat : number of different types of atoms in cell
2427 : !! typat(natom) : type of each atom
2428 : !! ngrid : number of gridpoints
2429 : !! r_vec_grid(3,ngrid) : real (non-reduced) coordinates for grid points
2430 : !! rho(ngrid) : input/output density array
2431 : !! a(3),b(3),c(3) : real-space basis vectors
2432 : !! atom_pos(3,natom) : reduced coordinates for atomic positions
2433 : !! natomgr(ntypat) : number of gridpoints for each atomic density grid
2434 : !! natomgrmax : max(natomgr(ntypat))
2435 : !! atomrgrid(natomgrmax,ntypat)
2436 : !! density(natomgrmax,ntypat)
2437 : !!
2438 : !! SIDE EFFECTS
2439 : !! rho(ngrid): input/output density array
2440 : !!
2441 : !! NOTES
2442 : !! There are two ways to compile the proto density in real space
2443 : !! for a solid. One alternative is that the density is calculated
2444 : !! for an extended grid encompassing the sphere of points around
2445 : !! one atom, and the results are folded back into the unit cell.
2446 : !! On the other hand one can, around each grid point, identify the
2447 : !! number of atoms in a sphere equivalent to the length of the radial
2448 : !! grid for each type of atom.
2449 : !! The second approach, with some modification, is taken here. The
2450 : !! numer of atoms in a supercell cell are listed such that the supercell
2451 : !! encompasses the atoms which could contribute to any point in the grid.
2452 : !! That list is kept and cycled through, to avoid recalculating it at
2453 : !! each point.
2454 : !! Note that the density calculated from the atom is the spherical
2455 : !! average, since there is no preferred direction without any
2456 : !! external field (and it's simpler)
2457 : !!
2458 : !! SOURCE
2459 :
2460 0 : subroutine atomden(MPI_enreg,natom,ntypat,typat,ngrid,r_vec_grid,rho,a,b,c,atom_pos, &
2461 0 : natomgr,natomgrmax,atomrgrid,density,prtvol,calctype)
2462 :
2463 : !Arguments ------------------------------------
2464 : !scalars
2465 : integer,intent(in) :: natom,ntypat,ngrid,natomgrmax,prtvol
2466 : character(len=7),intent(in) :: calctype
2467 : !arrays
2468 : type(MPI_type),intent(in) :: MPI_enreg
2469 : integer,intent(in) :: typat(natom),natomgr(ntypat)
2470 : real(dp),intent(in) :: r_vec_grid(3,ngrid),a(3),b(3),c(3)
2471 : real(dp),intent(in) :: atom_pos(3,natom),atomrgrid(natomgrmax,ntypat)
2472 : real(dp),intent(in) :: density(natomgrmax,ntypat)
2473 : real(dp),intent(inout) :: rho(ngrid)
2474 :
2475 : !Local variables-------------------------------
2476 : !scalars
2477 : character(len=500) :: message
2478 : integer :: cnt,delta,i,l,m,n,iatom,itypat,igrid,ncells,n_grid_p
2479 : integer :: ierr,spaceComm,nprocs,master,rank,remainder
2480 : real(dp) :: a_norm,b_norm,c_norm
2481 : real(dp) :: r_max,R_sphere_max,dp_dummy,ybcbeg,ybcend
2482 : !arrays
2483 0 : integer :: n_equiv_atoms(ntypat),grid_index(ngrid)
2484 0 : integer :: my_start_equiv_atoms(ntypat)
2485 0 : integer :: my_end_equiv_atoms(ntypat)
2486 0 : integer :: l_min(ntypat),m_min(ntypat),n_min(ntypat)
2487 0 : integer :: l_max(ntypat),m_max(ntypat),n_max(ntypat)
2488 : real(dp) :: center(3),dp_vec_dummy(3),delta_a(3),delta_b(3),delta_c(3)
2489 0 : real(dp) :: r_atom(3),grid_distances(ngrid)
2490 0 : integer, allocatable :: new_index(:),i_1d_dummy(:)
2491 0 : real(dp),allocatable :: equiv_atom_dist(:,:),equiv_atom_pos(:,:,:),rho_temp(:,:)
2492 0 : real(dp),allocatable :: dp_1d_dummy(:),dp_2d_dummy(:,:),ypp(:)
2493 0 : real(dp),allocatable :: x_fit(:),y_fit(:)
2494 : ! ************************************************************************
2495 :
2496 : !initialise and check parallel execution
2497 0 : spaceComm=MPI_enreg%comm_cell
2498 0 : nprocs=xmpi_comm_size(spaceComm)
2499 0 : rank=MPI_enreg%me_kpt
2500 :
2501 0 : master=0
2502 :
2503 : !initialise variables and vectors
2504 0 : a_norm = norm2(a)
2505 0 : b_norm = norm2(b)
2506 0 : c_norm = norm2(c)
2507 0 : center = (a+b+c)*half
2508 0 : dp_dummy = dot_product(a,b)/(b_norm*b_norm)
2509 0 : dp_vec_dummy = dp_dummy*b
2510 0 : delta_a = a - dp_vec_dummy
2511 0 : dp_dummy = dot_product(b,a)/(a_norm*a_norm)
2512 0 : dp_vec_dummy = dp_dummy*a
2513 0 : delta_b = b - dp_vec_dummy
2514 0 : dp_dummy = dot_product(c,(a+b))/(dot_product((a+b),(a+b)))
2515 0 : dp_vec_dummy = dp_dummy*(a+b)
2516 0 : delta_c = c - dp_vec_dummy
2517 0 : ABI_MALLOC(rho_temp,(ngrid,ntypat))
2518 0 : rho_temp = zero
2519 :
2520 : !write(std_out,*) '*** --- In atomden --- ***'
2521 : !write(std_out,*) ' a_norm:',a_norm,' b_norm:',b_norm,' c_norm:',c_norm
2522 : !write(std_out,*) 'delta_a:',delta_a,'delta_b:',delta_b,'delta_c:',delta_c
2523 : !write(std_out,*) ' center:',center
2524 :
2525 : !Find supercell which will contain all possible contributions
2526 : !for all atoms, and enumerate positions for all atoms
2527 : !TODO list of atoms can be "pruned", i.e identify all atoms
2528 : !that can't possibly contribute and remove from list.
2529 : !Should be most important for very oblique cells
2530 0 : do itypat=1,ntypat
2531 0 : R_sphere_max = atomrgrid(natomgr(itypat),itypat)
2532 0 : l_min(itypat) = -ceiling(R_sphere_max/norm2(delta_a))
2533 0 : l_max(itypat) = -l_min(itypat)
2534 0 : m_min(itypat) = -ceiling(R_sphere_max/norm2(delta_b))
2535 0 : m_max(itypat) = -m_min(itypat)
2536 0 : n_min(itypat) = -ceiling(R_sphere_max/norm2(delta_c))
2537 0 : n_max(itypat) = -n_min(itypat)
2538 : ncells = (l_max(itypat)-l_min(itypat)+1) &
2539 : & *(m_max(itypat)-m_min(itypat)+1) &
2540 0 : & *(n_max(itypat)-n_min(itypat)+1)
2541 0 : n_equiv_atoms(itypat) = 0
2542 0 : do iatom=1,natom
2543 0 : if (typat(iatom)==itypat) then
2544 0 : n_equiv_atoms(itypat) = n_equiv_atoms(itypat) + ncells
2545 : end if ! if type=itypat
2546 : end do ! number of atoms per cell
2547 0 : if ((rank==master).and.(prtvol>9)) then
2548 0 : write(message,'(a)') '*** --- In atomden --- find box ***'
2549 0 : call wrtout(std_out,message,'COLL')
2550 0 : write(message,'(a,I4)') ' itypat:',itypat
2551 0 : call wrtout(std_out,message,'COLL')
2552 0 : write(message,'(2(a,I4))') ' l_min:',l_min(itypat),' l_max:',l_max(itypat)
2553 0 : call wrtout(std_out,message,'COLL')
2554 0 : write(message,'(2(a,I4))') ' m_min:',m_min(itypat),' m_max:',m_max(itypat)
2555 0 : call wrtout(std_out,message,'COLL')
2556 0 : write(message,'(2(a,I4))') ' n_min:',n_min(itypat),' n_max:',n_max(itypat)
2557 0 : call wrtout(std_out,message,'COLL')
2558 0 : write(message,'(2(a,I4))') ' n_equiv_atoms:',n_equiv_atoms(itypat)
2559 0 : call wrtout(std_out,message,'COLL')
2560 : end if
2561 : end do !atom type
2562 :
2563 : !allocate arrays
2564 0 : n = maxval(n_equiv_atoms)
2565 0 : ABI_MALLOC(equiv_atom_pos,(3,n,ntypat))
2566 0 : ABI_MALLOC(equiv_atom_dist,(n,ntypat))
2567 0 : equiv_atom_pos = zero
2568 0 : equiv_atom_dist = zero
2569 :
2570 : !Find positions and distance of atoms from center of cell
2571 0 : do itypat=1,ntypat
2572 0 : i = 1
2573 0 : do l=l_min(itypat),l_max(itypat)
2574 0 : do m=m_min(itypat),m_max(itypat)
2575 0 : do n=n_min(itypat),n_max(itypat)
2576 0 : do iatom=1,natom
2577 0 : if (typat(iatom)==itypat) then
2578 0 : if (i>n_equiv_atoms(itypat)) then
2579 0 : ABI_ERROR('atomden: i>n_equiv_atoms')
2580 : end if
2581 : equiv_atom_pos(:,i,itypat) = (atom_pos(1,iatom)+dble(l))*a &
2582 : & + (atom_pos(2,iatom)+dble(m))*b &
2583 0 : & + (atom_pos(3,iatom)+dble(n))*c
2584 0 : dp_vec_dummy = equiv_atom_pos(:,i,itypat)-center
2585 0 : equiv_atom_dist(i,itypat) = norm2(dp_vec_dummy)
2586 0 : i = i + 1
2587 : end if
2588 : end do
2589 : end do !n
2590 : end do !m
2591 : end do !l
2592 : ! write(std_out,*) '*** --- In atomden --- find equiv ***'
2593 : ! write(std_out,*) ' itypat:',itypat
2594 : ! write(std_out,*) ' equiv_atom_pos:'
2595 : ! write(std_out,*) equiv_atom_pos(:,:,itypat)
2596 : ! write(std_out,*) ' equiv_atom_dist:',equiv_atom_dist(:,itypat)
2597 : end do !atom type
2598 :
2599 : !Sort the atoms after distance so that the density from the ones
2600 : !furthest away can be added first. This is to prevent truncation error.
2601 0 : do itypat=1,ntypat
2602 0 : n = n_equiv_atoms(itypat)
2603 0 : ABI_MALLOC(dp_1d_dummy,(n))
2604 0 : ABI_MALLOC(new_index,(n))
2605 0 : ABI_MALLOC(dp_2d_dummy,(3,n))
2606 0 : dp_1d_dummy = equiv_atom_dist(1:n,itypat)
2607 0 : dp_2d_dummy = equiv_atom_pos(1:3,1:n,itypat)
2608 0 : do i=1,n
2609 0 : new_index(i) = i
2610 : end do
2611 0 : call sort_dp(n,dp_1d_dummy,new_index,tol14)
2612 0 : do i=1,n
2613 : !write(std_out,*) i,' -> ',new_index(i)
2614 0 : equiv_atom_pos(1:3,n+1-i,itypat) = dp_2d_dummy(1:3,new_index(i))
2615 0 : equiv_atom_dist(1:n,itypat) = dp_1d_dummy
2616 : end do
2617 0 : ABI_FREE(dp_1d_dummy)
2618 0 : ABI_FREE(new_index)
2619 0 : ABI_FREE(dp_2d_dummy)
2620 : ! write(std_out,*) '*** --- In atomden --- sorting atoms ***'
2621 : ! write(std_out,*) ' itypat:',itypat
2622 : ! write(std_out,*) ' equiv_atom_pos:'
2623 : ! write(std_out,*) equiv_atom_pos(:,:,itypat)
2624 : ! write(std_out,*) ' equiv_atom_dist:',equiv_atom_dist(:,itypat)
2625 : end do ! atom type
2626 :
2627 : !Divide the work in case of parallel execution
2628 0 : if (nprocs==1) then ! Make sure everything runs with one proc
2629 0 : if (prtvol>9) then
2630 0 : write(message,'(a)') ' In atomden - number of processors: 1'
2631 0 : call wrtout(std_out,message,'COLL')
2632 0 : write(message,'(a)') ' Calculation of proto-atomic density done in serial'
2633 0 : call wrtout(std_out,message,'COLL')
2634 : end if
2635 0 : do itypat=1,ntypat
2636 0 : if (prtvol>9) then
2637 0 : write(message,'(a,I6)') ' Number of equivalent atoms:',n_equiv_atoms(itypat)
2638 0 : call wrtout(std_out,message,'COLL')
2639 : end if
2640 0 : my_start_equiv_atoms(itypat) = 1
2641 0 : my_end_equiv_atoms(itypat) = n_equiv_atoms(itypat)
2642 : end do
2643 : else
2644 0 : if (rank==master.and.prtvol>9) then
2645 0 : write(message,'(a,I5)') ' In atomden - number of processors:',nprocs
2646 0 : call wrtout(std_out,message,'COLL')
2647 0 : write(message,'(a)') ' Calculation of proto-atomic density done in parallel'
2648 0 : call wrtout(std_out,message,'COLL')
2649 : end if
2650 0 : do itypat=1,ntypat
2651 0 : if (rank==master.and.prtvol>9) then
2652 0 : write(message,'(a,I6)') ' Number of equivalent atoms:',n_equiv_atoms(itypat)
2653 0 : call wrtout(std_out,message,'COLL')
2654 : end if
2655 : ! Divide the atoms among the processors by shuffling indices
2656 0 : delta = int(floor(real(n_equiv_atoms(itypat))/real(nprocs)))
2657 0 : remainder = n_equiv_atoms(itypat)-nprocs*delta
2658 0 : my_start_equiv_atoms(itypat) = 1+rank*delta
2659 0 : my_end_equiv_atoms(itypat) = (rank+1)*delta
2660 : ! Divide the remainder points among the processors
2661 : ! by shuffling indices
2662 0 : if ((rank+1)>remainder) then
2663 0 : my_start_equiv_atoms(itypat) = my_start_equiv_atoms(itypat) + remainder
2664 0 : my_end_equiv_atoms(itypat) = my_end_equiv_atoms(itypat) + remainder
2665 : else
2666 0 : my_start_equiv_atoms(itypat) = my_start_equiv_atoms(itypat) + rank
2667 0 : my_end_equiv_atoms(itypat) = my_end_equiv_atoms(itypat) + rank + 1
2668 : end if
2669 0 : if (prtvol>9) then
2670 0 : write(message,'(a,I3)') ' For atom type: ',itypat
2671 0 : call wrtout(std_out,message,'PERS')
2672 : ! write(message,'(a,I6)') ' I''ll take atoms from: ',my_start_equiv_atoms(itypat)
2673 : ! call wrtout(std_out,message,'PERS')
2674 : ! write(message,'(a,I6)') ' total for me: ',my_end_equiv_atoms(itypat)
2675 : ! call wrtout(std_out,message,'PERS')
2676 0 : write(message,'(a,I6)') ' total for me: ', &
2677 0 : & my_end_equiv_atoms(itypat)+1-my_start_equiv_atoms(itypat)
2678 0 : call wrtout(std_out,message,'PERS')
2679 : end if
2680 : end do
2681 : end if
2682 :
2683 : !Loop over types of atoms and equivalent atoms and
2684 : !interpolate density onto grid
2685 0 : do itypat=1,ntypat
2686 : ! do iatom=my_start_equiv_atoms(itypat),my_end_equiv_atoms(itypat)
2687 :
2688 0 : cnt = 0
2689 0 : iatom = rank+1 - nprocs
2690 : ! Parallel execution of loop
2691 0 : do
2692 0 : cnt = cnt + 1
2693 0 : iatom = iatom + nprocs
2694 0 : if (iatom>n_equiv_atoms(itypat)) exit ! Exit if index is too large
2695 :
2696 0 : if (mod(cnt,100)==0.and.prtvol>0) then
2697 0 : write(message,'(2(a,I6))') ' atoms so far',cnt,' of: ',n_equiv_atoms(itypat)/nprocs
2698 0 : call wrtout(std_out,message,'PERS')
2699 : end if
2700 :
2701 0 : r_max = atomrgrid(natomgr(itypat),itypat)
2702 0 : r_atom = equiv_atom_pos(:,iatom,itypat)
2703 :
2704 : ! Set up an array with the gridpoint distances
2705 0 : i = 1
2706 0 : grid_distances = zero
2707 0 : grid_index = 0
2708 0 : do igrid=1,ngrid
2709 0 : dp_vec_dummy(:) = r_vec_grid(:,igrid) - r_atom(:)
2710 0 : dp_dummy = norm2(dp_vec_dummy)
2711 0 : if (dp_dummy <= r_max) then
2712 0 : grid_distances(i) = dp_dummy
2713 0 : grid_index(i) = igrid
2714 0 : i = i + 1
2715 : else
2716 0 : cycle ! cycle if point is too far away
2717 : end if
2718 : end do
2719 0 : n_grid_p = i - 1
2720 :
2721 0 : if (n_grid_p==0) cycle ! Cycle if no point needs
2722 : ! to be interpolated
2723 :
2724 : ! Sort points to be interpolated in ascending order
2725 0 : ABI_MALLOC(dp_1d_dummy,(n_grid_p))
2726 0 : ABI_MALLOC(new_index,(n_grid_p))
2727 0 : ABI_MALLOC(i_1d_dummy,(n_grid_p))
2728 0 : dp_1d_dummy = grid_distances(1:n_grid_p)
2729 0 : do i=1,n_grid_p
2730 0 : new_index(i) = i
2731 : end do
2732 0 : call sort_dp(n_grid_p,dp_1d_dummy,new_index,tol16)
2733 0 : grid_distances(1:n_grid_p) = dp_1d_dummy
2734 0 : i_1d_dummy = grid_index(1:n_grid_p)
2735 0 : do i=1,n_grid_p
2736 : !write(std_out,*) i_1d_dummy(i),' -> ',i_1d_dummy(new_index(i))
2737 0 : grid_index(i) = i_1d_dummy(new_index(i))
2738 : end do
2739 0 : ABI_FREE(dp_1d_dummy)
2740 0 : ABI_FREE(new_index)
2741 0 : ABI_FREE(i_1d_dummy)
2742 :
2743 : ! Interpolate density onto all grid points
2744 0 : ABI_MALLOC(ypp,(natomgr(itypat)))
2745 0 : ABI_MALLOC(x_fit,(n_grid_p))
2746 0 : ABI_MALLOC(y_fit,(n_grid_p))
2747 0 : ypp = zero; y_fit = zero
2748 0 : ybcbeg = zero; ybcend = zero
2749 0 : x_fit = grid_distances(1:n_grid_p)
2750 : call spline(atomrgrid(1:natomgr(itypat),itypat), &
2751 : & density(1:natomgr(itypat),itypat), &
2752 0 : & natomgr(itypat),ybcbeg,ybcend,ypp)
2753 : call splint(natomgr(itypat),atomrgrid(1:natomgr(itypat),itypat), &
2754 : & density(1:natomgr(itypat),itypat),ypp,n_grid_p, &
2755 0 : & x_fit,y_fit)
2756 :
2757 : ! Save the interpolated points to grid
2758 0 : do i=1,n_grid_p
2759 0 : rho_temp(grid_index(i),itypat) = rho_temp(grid_index(i),itypat) + y_fit(i)
2760 : end do
2761 0 : ABI_FREE(ypp)
2762 0 : ABI_FREE(x_fit)
2763 0 : ABI_FREE(y_fit)
2764 :
2765 : end do ! n equiv atoms
2766 : end do ! type of atom
2767 :
2768 : ! Collect all contributions to rho_temp if we are running in parallel
2769 0 : if (nprocs>1) then
2770 0 : call xmpi_barrier(spaceComm)
2771 0 : call xmpi_sum_master(rho_temp,master,spaceComm,ierr)
2772 0 : call xmpi_barrier(spaceComm)
2773 0 : if (prtvol>9) then
2774 0 : write(message,'(a)') ' In atomden - contributions to rho_temp collected'
2775 0 : call wrtout(std_out,message,'PERS')
2776 : end if
2777 : end if
2778 :
2779 : !Now rho_temp contains the atomic protodensity for each atom.
2780 : !Check whether this is to replace or be added to the input/output array
2781 : !and sum up contributions
2782 0 : if (trim(calctype)=='replace') rho = zero
2783 0 : do itypat=1,ntypat
2784 0 : rho(:) = rho(:) + rho_temp(:,itypat)
2785 : end do
2786 :
2787 : ! deallocations
2788 0 : ABI_SFREE(rho_temp)
2789 0 : ABI_SFREE(equiv_atom_pos)
2790 0 : ABI_SFREE(equiv_atom_dist)
2791 :
2792 0 : end subroutine atomden
2793 : !!***
2794 :
2795 : end module m_mkrho
2796 : !!***
|
