Branch data 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 : 17609 : subroutine mkrho(cg,dtset,gprimd,irrzon,kg,mcg,mpi_enreg,npwarr,occ,paw_dmft,phnons,&
139 : 17609 : & 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 : 17609 : 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 : 17609 : real(dp),allocatable :: cwavef_rot(:,:,:,:),occ_diag(:),occ_k(:)
188 : 17609 : 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 : 17609 : call timab(790+tim_mkrho,1,tsec)
208 : 17609 : call timab(799,1,tsec)
209 : :
210 [ + + ]: 17609 : if(mpi_enreg%paralbd==0) tim_fourwf=3
211 [ + + ]: 17609 : if(mpi_enreg%paralbd==1) tim_fourwf=6
212 : :
213 [ + + ]: 17609 : if(.not.(present(option))) then
214 : 16767 : ioption=0
215 : : else
216 : 842 : ioption=option
217 : : end if
218 : :
219 [ + + ]: 17609 : if(.not.(present(printout))) then
220 : : l_printout = .true.
221 : : else
222 : 1593 : l_printout = printout
223 : : end if
224 : :
225 : : ! Not sure what to do for Wannier90 DMFT
226 [ + + - + ]: 17609 : 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 [ + + ]: 17609 : 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 [ - + ]: 17471 : else if(paw_dmft%use_sc_dmft==10) then
232 : 0 : nbandc1=paw_dmft%mbandc
233 : : else
234 : : nbandc1=1
235 : : end if
236 : 17609 : 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 : 17609 : my_nspinor=max(1,dtset%nspinor/mpi_enreg%nproc_spinor)
246 [ + + ]: 17609 : if (mpi_enreg%paral_spinor==0) then
247 : 17605 : ispinor_index=1;jspinor_index=1
248 : 17605 : nspinor1TreatedByThisProc=.true.
249 : 17605 : 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 [ + - - + ]: 17609 : ABI_BUG(sjoin('ioption argument value should be 0,1 or 2 while got:', itoa(ioption)))
275 : : end select
276 : :
277 : : !Init me
278 : 17609 : me=mpi_enreg%me_kpt
279 : : !zero the charge density array in real space
280 : : !$OMP PARALLEL DO COLLAPSE(2)
281 [ + + ]: 43170 : do ispden=1,dtset%nspden
282 [ + + ]: 241581654 : do ifft=1,dtset%nfft
283 : 241564045 : rhor(ifft,ispden)=zero
284 : : end do
285 : : end do
286 : :
287 : : !WVL - Branching with a separate mkrho procedure in wavelet.
288 [ - + ]: 17609 : 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 : 17609 : gpu_option=dtset%gpu_option
304 [ - + - - ]: 17609 : gpu_cwavef=(gpu_option==ABI_GPU_OPENMP .and. paw_dmft%use_sc_dmft/=1)
305 : 17609 : 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 [ - + ]: 17609 : 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 [ + + ]: 36866 : do alpha=1,nalpha
321 [ + + ]: 56123 : do beta=1,nbeta
322 : :
323 : : ! start loop over spin and k points
324 : 19257 : bdtot_index=0
325 : 19257 : icg=0
326 : :
327 : : ! n4,n5,n6 are FFT dimensions, modified to avoid cache trashing
328 : 19257 : n1 = dtset%ngfft(1) ; n2 = dtset%ngfft(2) ; n3 = dtset%ngfft(3)
329 : 19257 : n4 = dtset%ngfft(4) ; n5 = dtset%ngfft(5) ; n6 = dtset%ngfft(6)
330 : :
331 : 19257 : ndat = 1
332 [ + + ]: 19257 : if (mpi_enreg%paral_kgb==1) then
333 : 4488 : ndat = mpi_enreg%bandpp
334 [ - + ]: 14769 : else if (gpu_option/=ABI_GPU_DISABLED) then
335 : 0 : ndat = dtset%mband
336 : : end if
337 : :
338 [ + - ]: 19257 : 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 [ - + ]: 19257 : if (gpu_option/=ABI_GPU_DISABLED) then
344 [ # # # # : 0 : ABI_MALLOC(cwavef,(2,dtset%mpw*ndat,my_nspinor))
# # # # #
# # # # #
# # ]
345 : : else
346 [ - + + - : 77028 : ABI_MALLOC(cwavef,(2,dtset%mpw,my_nspinor))
- + + - +
- + - - +
- + ]
347 : : end if
348 : : end if
349 : :
350 [ + - ]: 19257 : 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 [ - + + - : 96285 : 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 [ - + + - : 96285 : ABI_MALLOC(wfraug, (2,n4,n5,n6*(ndat/l_nfft_blocks+ndat-(ndat/l_nfft_blocks)*l_nfft_blocks)))
- + + - -
+ + - + -
+ - - + -
+ ]
372 [ + + + - : 77028 : ABI_MALLOC(cwavefb, (2,dtset%mpw*paw_dmft%use_sc_dmft,my_nspinor))
- + + - +
- + + - +
- + ]
373 [ + + ]: 19257 : if(dtset%nspden==4) then
374 [ - + + - : 8092 : ABI_MALLOC(rhoaug_up, (n4,n5,n6))
- + + - -
+ + - - +
- + ]
375 [ - + + - : 8092 : ABI_MALLOC(rhoaug_down,(n4,n5,n6))
- + + - -
+ + - - +
- + ]
376 [ - + + - : 8092 : ABI_MALLOC(rhoaug_mx, (n4,n5,n6))
- + + - -
+ + - - +
- + ]
377 [ - + + - : 8092 : ABI_MALLOC(rhoaug_my, (n4,n5,n6))
- + + - -
+ + - - +
- + ]
378 [ + + + + : 15871707 : rhoaug_up(:,:,:)=zero
+ + ]
379 [ + + + + : 15871707 : rhoaug_down(:,:,:)=zero
+ + ]
380 [ + + + + : 15871707 : rhoaug_mx(:,:,:)=zero
+ + ]
381 [ + + + + : 15871707 : rhoaug_my(:,:,:)=zero
+ + ]
382 : : end if
383 : : end if
384 : :
385 [ + + ]: 40698 : do isppol=1,dtset%nsppol
386 : 21441 : ikg=0
387 : :
388 [ - + ]: 21441 : 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 [ + + + + : 311170954 : rhoaug(:,:,:)=zero
+ + ]
392 : : end if
393 [ + + ]: 252718 : do ikpt=1,dtset%nkpt
394 : :
395 : 231277 : nband_k = dtset%nband(ikpt+(isppol-1)*dtset%nkpt)
396 : 231277 : mband_mem = nband_k
397 [ + + ]: 231277 : if (dtset%paral_kgb==0) mband_mem = nband_k/mpi_enreg%nproc_band
398 : 231277 : npw_k=npwarr(ikpt)
399 : 231277 : istwf_k = dtset%istwfk(ikpt)
400 : :
401 [ + + ]: 231277 : 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 [ + - - + : 740424 : ABI_MALLOC(gbound,(2*dtset%mgfft+8,2))
+ - + - +
- + - - +
- + - + ]
407 [ + - ]: 185106 : 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 [ - + + - : 555318 : ABI_MALLOC(kg_k,(3,npw_k))
+ - + - -
+ - + ]
413 : : end if
414 : :
415 [ + + + + ]: 158285526 : kg_k(:,1:npw_k)=kg(:,1+ikg:npw_k+ikg)
416 : 185106 : 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 [ + + ]: 185106 : 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 [ - + - - ]: 169098 : 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 : 169098 : iband_me = 0
556 [ + + ]: 1760759 : do iband=1,nband_k
557 [ + + ]: 1591661 : if(mpi_enreg%paralbd==1)then
558 [ + + ]: 259639 : if(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,iband,iband,isppol,me)) cycle
559 : : end if
560 : 1591505 : iband_me = iband_me + 1
561 [ + + ]: 3543832 : do ibandc1=1,nbandc1 ! in case of DMFT
562 : : ! Check if DMFT and only treat occupied states (check on occ.)
563 [ + + ]: 1783229 : 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 : 1576279 : use_nondiag_occup_dmft = 0
577 : 1576279 : locc_test = abs(occ(iband+bdtot_index))>tol8
578 : : end if
579 [ + + ]: 3339546 : 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 : 1016807 : ipwsp=(iband_me-1)*npw_k*my_nspinor +icg
584 [ + - + + : 1682272277 : cwavef(:,1:npw_k,1) = cg(:,1+ipwsp :ipwsp+npw_k)
+ + + + +
+ ]
585 [ + + + + : 256826597 : if (my_nspinor==2) cwavef(:,1:npw_k,2)=cg(:,1+ipwsp+npw_k:ipwsp+2*npw_k)
+ + + + +
+ ]
586 [ + + ]: 1016807 : 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 [ - + ]: 990230 : 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 [ + + ]: 1016807 : 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 : 988699 : weight=occ(iband+bdtot_index)*dtset%wtk(ikpt)/ucvol
611 : 988699 : 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 : 1016807 : & gpu_option=gpu_option)
621 [ + + ]: 1016807 : if(dtset%nspinor==2)then
622 [ + + ]: 181273 : 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 [ - + ]: 140109 : 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 [ + - ]: 140109 : 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 [ - + + - : 420327 : ABI_MALLOC(cwavef_x,(2,npw_k))
+ - + - -
+ - + ]
643 [ - + + - : 280218 : ABI_MALLOC(cwavef_y,(2,npw_k))
- + - + ]
644 [ + + + - : 420327 : ABI_MALLOC(cwavefb_x,(2,npw_k*paw_dmft%use_sc_dmft))
+ - + + -
+ - + ]
645 [ + + + - : 280218 : ABI_MALLOC(cwavefb_y,(2,npw_k*paw_dmft%use_sc_dmft))
- + - + ]
646 : : end if
647 : : ! $(\Psi^{1}+\Psi^{2})$
648 [ + - + + : 189368043 : cwavef_x(:,:)=cwavef(:,1:npw_k,1)+cwavef(:,1:npw_k,2)
+ + + + +
+ ]
649 : : ! $(\Psi^{1}-i \Psi^{2})$
650 [ + - + + : 63216087 : cwavef_y(1,:)=cwavef(1,1:npw_k,1)+cwavef(2,1:npw_k,2)
+ + ]
651 [ + + + + ]: 63216087 : cwavef_y(2,:)=cwavef(2,1:npw_k,1)-cwavef(1,1:npw_k,2)
652 [ + + ]: 140109 : 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 [ + - + + : 1820319987 : 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 : 140109 : & 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 : 140109 : & 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 : 140109 : & gpu_option=gpu_option)
678 [ - + ]: 140109 : 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 [ - + ]: 140109 : ABI_FREE(cwavef_x)
687 [ - + ]: 140109 : ABI_FREE(cwavef_y)
688 [ - + ]: 140109 : ABI_FREE(cwavefb_x)
689 [ - + ]: 140109 : 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 : 731078 : 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 [ - + ]: 185106 : ABI_FREE(gbound)
979 : :
980 [ + - ]: 185106 : 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 [ - + ]: 185106 : ABI_FREE(kg_k)
986 : : end if
987 : :
988 : 185106 : bdtot_index=bdtot_index+nband_k
989 : :
990 [ + - ]: 206547 : if (dtset%mkmem/=0) then
991 : 185106 : icg=icg+npw_k*my_nspinor*mband_mem !iband_me
992 : 185106 : 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 [ + + ]: 21441 : 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 [ + + ]: 40698 : if(dtset%nspden==1 .or. dtset%nspden==2) then
1017 : 19418 : call fftpac(isppol,mpi_enreg,dtset%nspden,n1,n2,n3,n4,n5,n6,dtset%ngfft,rhor,rhoaug,1)
1018 [ + - ]: 2023 : else if(dtset%nspden==4) then
1019 : 2023 : ispden=1
1020 : 2023 : call fftpac(ispden,mpi_enreg,dtset%nspden,n1,n2,n3,n4,n5,n6,dtset%ngfft,rhor,rhoaug_up,1)
1021 : 2023 : ispden=2
1022 : 2023 : call fftpac(ispden,mpi_enreg,dtset%nspden,n1,n2,n3,n4,n5,n6,dtset%ngfft,rhor,rhoaug_mx,1)
1023 : 2023 : ispden=3
1024 : 2023 : call fftpac(ispden,mpi_enreg,dtset%nspden,n1,n2,n3,n4,n5,n6,dtset%ngfft,rhor,rhoaug_my,1)
1025 : 2023 : ispden=4
1026 : 2023 : call fftpac(ispden,mpi_enreg,dtset%nspden,n1,n2,n3,n4,n5,n6,dtset%ngfft,rhor,rhoaug_down,1)
1027 [ + - ]: 2023 : 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 [ - + ]: 2023 : ABI_FREE(rhoaug_up)
1036 [ - + ]: 2023 : ABI_FREE(rhoaug_down)
1037 [ - + ]: 2023 : ABI_FREE(rhoaug_mx)
1038 [ - + ]: 2023 : ABI_FREE(rhoaug_my)
1039 : : end if
1040 : : end if
1041 : :
1042 : : end do ! isppol
1043 : :
1044 [ + - ]: 19257 : 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 [ + + ]: 19257 : if(associated(cwavef)) then
1057 : 14769 : ABI_FREE(cwavef)
1058 : : end if
1059 [ + - ]: 19257 : if(associated(cwavefb)) then
1060 : 19257 : 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 [ - + ]: 19257 : ABI_FREE(rhoaug)
1066 [ - + ]: 19257 : ABI_FREE(wfraug)
1067 : : end if
1068 : :
1069 [ + + ]: 19257 : 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 : 19257 : call timab(48,1,tsec)
1077 : 19257 : call timab(71,1,tsec)
1078 : 19257 : spaceComm=mpi_enreg%comm_cell
1079 [ - + ]: 19257 : if (mpi_enreg%paral_hf==1)spaceComm=mpi_enreg%comm_kpt
1080 [ + + ]: 19257 : if(mpi_enreg%paral_kgb==1)spaceComm=mpi_enreg%comm_kpt
1081 : 19257 : call xmpi_sum(rhor,spaceComm,ierr)
1082 : 19257 : call timab(71,2,tsec)
1083 : 19257 : call timab(48,2,tsec)
1084 : :
1085 : 19257 : call timab(799,2,tsec)
1086 : 19257 : call timab(549,1,tsec)
1087 : :
1088 [ + + ]: 38514 : 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 [ + + ]: 17609 : 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 : 17609 : 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 [ + + ]: 17609 : if(present(extfpmd)) then
1125 [ - + ]: 8361 : 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 : 17609 : 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 : 17609 : phnons,rhog,rhor,rprimd,dtset%symafm,dtset%symrel,dtset%tnons)
1138 [ + + ]: 17609 : 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 [ + - - ]: 17609 : 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 : 17609 : 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 : 17609 : call timab(799,1,tsec)
1160 : :
1161 [ + + ]: 17609 : 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 [ + + ]: 17609 : if (l_printout) then
1168 : 16016 : call wrtout(std_out,' mkrho: echo density (plane-wave part only)','COLL')
1169 : 16016 : call prtrhomxmn(std_out,mpi_enreg,dtset%nfft,dtset%ngfft,dtset%nspden,1,rhor,optrhor=ioption,ucvol=ucvol)
1170 : : end if
1171 : :
1172 : 17609 : call timab(799,2,tsec)
1173 : 17609 : call timab(790+tim_mkrho,2,tsec)
1174 : :
1175 : : DBG_EXIT("COLL")
1176 : :
1177 [ - + - + : 17609 : 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 : 3096 : subroutine initro(atindx,densty,gmet,gsqcut,izero,mgfft,mpi_enreg,mqgrid,natom,nattyp,&
1222 : 3096 : & 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 : 3096 : integer, contiguous, pointer :: fftn2_distrib(:),ffti2_local(:),fftn3_distrib(:),ffti3_local(:)
1251 : 6192 : real(dp) :: length(ntypat)
1252 : 3096 : real(dp),allocatable :: work(:), spinat_indx(:,:)
1253 : 6192 : logical :: use_gaussian(ntypat)
1254 : : ! *************************************************************************
1255 : :
1256 : 3096 : if (nspden==4) then
1257 : 124 : 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 : 3096 : n1=ngfft(1); n2=ngfft(2); n3=ngfft(3); me_fft=ngfft(11); nproc_fft=ngfft(10)
1262 : :
1263 [ - + + - : 9288 : ABI_MALLOC(work,(nfft))
+ - + - -
+ - + ]
1264 [ - + + - : 9288 : ABI_MALLOC(spinat_indx,(3,natom))
+ - + - -
+ - + ]
1265 : :
1266 : : ! Get the distrib associated with this fft_grid
1267 : 3096 : 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 [ + + ]: 10850 : do ia=1,natom
1271 [ + + ]: 34112 : spinat_indx(:,atindx(ia))=spinat(:,ia)
1272 : : end do
1273 : :
1274 : : ! Check whether the values of spinat are acceptable
1275 [ + + ]: 3096 : if (nspden==2)then
1276 : 386 : ia1=1
1277 [ + + ]: 831 : do itypat=1,ntypat
1278 : : !ia1,ia2 sets range of loop over atoms:
1279 : 445 : ia2=ia1+nattyp(itypat)-1
1280 [ + + ]: 1088 : do ia=ia1,ia2
1281 [ - + ]: 643 : if( sqrt(spinat_indx(1,ia)**2+spinat_indx(2,ia)**2+spinat_indx(3,ia)**2) &
1282 : 445 : > 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 : 831 : 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 : 3096 : jtemp=0
1296 [ + + ]: 7480 : do itypat=1,ntypat
1297 : 4384 : use_gaussian(itypat)=.true.
1298 [ + + ]: 4384 : if (usepaw==0) use_gaussian(itypat) = .not. psps%nctab(itypat)%has_tvale
1299 [ + + ]: 4384 : if (usepaw==1) use_gaussian(itypat)=(pawtab(itypat)%has_tvale==0)
1300 [ + + ]: 4384 : if (.not.use_gaussian(itypat)) jtemp=jtemp+1
1301 : :
1302 [ + + ]: 7480 : if (use_gaussian(itypat)) then
1303 : 2920 : 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 : 2920 : ' initro: for itypat=',itypat,', take decay length=',length(itypat),',',ch10,&
1306 : 5840 : ' initro: indeed, coreel=',znucl(itypat)-zion(itypat),', nval=',int(zion(itypat)),' and densty=',densty(itypat,1),'.'
1307 : 2920 : call wrtout(std_out,message)
1308 : : else
1309 : 1464 : write(message,"(a,i3,a)")' initro: for itypat=',itypat,", take pseudo charge density from pp file"
1310 : 1464 : call wrtout(std_out,message)
1311 : : end if
1312 : : end do
1313 : :
1314 [ + + ]: 3096 : if (jtemp>0) then
1315 : 940 : dq=(qgrid(mqgrid)-qgrid(1))/dble(mqgrid-1)
1316 : 940 : dqm1=1.0_dp/dq
1317 : 940 : dq2div6=dq**2/6.0_dp
1318 : : end if
1319 : :
1320 : 3096 : cutoff=gsqcut*tolfix
1321 : 3096 : xnorm=1.0_dp/ucvol
1322 : :
1323 : 3096 : id1=n1/2+2
1324 : 3096 : id2=n2/2+2
1325 : 3096 : id3=n3/2+2
1326 : :
1327 [ + + ]: 3096 : if(nspden /= 4) then
1328 : :
1329 [ + + ]: 6330 : 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 [ + + + + ]: 294658825 : rhog = zero
1337 : :
1338 : 3358 : ia1=1
1339 [ + + ]: 8031 : do itypat=1,ntypat
1340 [ + + ]: 4673 : if (use_gaussian(itypat)) alf2pi2=(two_pi*length(itypat))**2
1341 : :
1342 : : ! ia1,ia2 sets range of loop over atoms:
1343 : 4673 : 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 [ + + ]: 118220 : do i3=1,n3
1348 : 113547 : ig3=i3-(i3/id3)*n3-1
1349 [ + + ]: 3451291 : do i2=1,n2
1350 : 3333071 : ig2=i2-(i2/id2)*n2-1
1351 [ + + ]: 3446618 : if (fftn2_distrib(i2)==me_fft) then
1352 [ + + ]: 143999561 : do i1=1,n1
1353 : :
1354 : 140807298 : ig1=i1-(i1/id1)*n1-1
1355 : 140807298 : 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 : 140807298 : dble(2*ig2*ig3)*gmet(2,3)+dble(2*ig3*ig1)*gmet(3,1)
1359 : :
1360 : : ! Skip G**2 outside cutoff:
1361 [ + + ]: 143999561 : 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 : 57506938 : sfr=zero;sfi=zero
1367 [ + + ]: 57506938 : if (ispden==1) then
1368 [ + + ]: 151333671 : do ia=ia1,ia2
1369 : 104715927 : sfr=sfr+phre_ini(ig1,ig2,ig3,ia)
1370 : 151333671 : sfi=sfi-phimag_ini(ig1,ig2,ig3,ia)
1371 : : end do
1372 [ + + ]: 46617744 : if (use_gaussian(itypat)) then
1373 : 25178408 : sfr=sfr*zion(itypat)
1374 : 25178408 : sfi=sfi*zion(itypat)
1375 : : end if
1376 : : else
1377 [ + + ]: 10889194 : fact0=half;if (.not.use_gaussian(itypat)) fact0=half/zion(itypat)
1378 [ + + ]: 21763934 : do ia=ia1,ia2
1379 : : ! Here, take care only of the z component
1380 : 10874740 : fact=fact0*(zion(itypat)+spinat_indx(3,ia))
1381 : 10874740 : sfr=sfr+phre_ini(ig1,ig2,ig3,ia)*fact
1382 : 21763934 : sfi=sfi-phimag_ini(ig1,ig2,ig3,ia)*fact
1383 : : end do
1384 : : end if
1385 : :
1386 : : ! Charge density integrating to one
1387 [ + + ]: 57506938 : if (use_gaussian(itypat)) then
1388 : 34478884 : rhoat=xnorm*exp(-gsquar*alf2pi2)
1389 : : ! Multiply structure factor times rhoat (atomic density in reciprocal space)
1390 : 34478884 : rhog(re,ii)=rhog(re,ii)+sfr*rhoat
1391 : 34478884 : rhog(im,ii)=rhog(im,ii)+sfi*rhoat
1392 : : else
1393 : 23028054 : gmag=sqrt(gsquar)
1394 : 23028054 : jj=1+int(gmag*dqm1)
1395 : 23028054 : diff=gmag-qgrid(jj)
1396 : 23028054 : bb = diff*dqm1
1397 : 23028054 : aa = one-bb
1398 : 23028054 : cc = aa*(aa**2-one)*dq2div6
1399 : 23028054 : dd = bb*(bb**2-one)*dq2div6
1400 [ + + ]: 23028054 : if (usepaw == 1) then
1401 : : rhoat=(aa*pawtab(itypat)%tvalespl(jj,1)+bb*pawtab(itypat)%tvalespl(jj+1,1)+&
1402 : 17771203 : & cc*pawtab(itypat)%tvalespl(jj,2)+dd*pawtab(itypat)%tvalespl(jj+1,2)) *xnorm
1403 [ + - ]: 5256851 : else if (usepaw == 0) then
1404 : : rhoat=(aa*psps%nctab(itypat)%tvalespl(jj,1)+bb*psps%nctab(itypat)%tvalespl(jj+1,1)+&
1405 : 5256851 : 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 : 23028054 : rhog(re,ii)=rhog(re,ii)+sfr*rhoat
1411 : 23028054 : 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 : 8031 : ia1=ia2+1
1421 : :
1422 : : end do ! itypat
1423 : :
1424 : : ! Set contribution of unbalanced components to zero
1425 [ + + ]: 3358 : if (izero==1) then
1426 : 975 : 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 : 3358 : call fourdp(1,rhog,work,1,mpi_enreg,nfft,1,ngfft,0)
1432 [ + + ]: 98224819 : rhor(:,ispden)=work(:)
1433 : : end do ! End loop on spins
1434 : :
1435 : : else if (nspden==4) then
1436 : :
1437 [ + + ]: 620 : 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 [ + + + + ]: 19919248 : rhog(:,:)=zero
1445 : :
1446 : 496 : ia1=1
1447 [ + + ]: 1120 : do itypat=1,ntypat
1448 : :
1449 [ + + ]: 624 : if (use_gaussian(itypat)) alf2pi2=(two_pi*length(itypat))**2
1450 : :
1451 : : ! ia1,ia2 sets range of loop over atoms:
1452 : 624 : ia2=ia1+nattyp(itypat)-1
1453 : 624 : ii=0
1454 [ + + ]: 15060 : do i3=1,n3
1455 : 14436 : ig3=i3-(i3/id3)*n3-1
1456 [ + + ]: 347472 : do i2=1,n2
1457 : 332412 : ig2=i2-(i2/id2)*n2-1
1458 [ + + ]: 346848 : if (fftn2_distrib(i2)==me_fft) then
1459 [ + + ]: 8776104 : do i1=1,n1
1460 : :
1461 : 8452908 : ig1=i1-(i1/id1)*n1-1
1462 : 8452908 : 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 : 8452908 : dble(2*ig2*ig3)*gmet(2,3)+dble(2*ig3*ig1)*gmet(3,1)
1466 : :
1467 : : ! Skip G**2 outside cutoff:
1468 [ + + ]: 8776104 : 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 : 2831520 : sfr=zero;sfi=zero
1473 [ + + ]: 2831520 : if(ispden==1)then
1474 [ + + ]: 2030007 : do ia=ia1,ia2
1475 : 1322127 : sfr=sfr+phre_ini(ig1,ig2,ig3,ia)
1476 : 2030007 : sfi=sfi-phimag_ini(ig1,ig2,ig3,ia)
1477 : : end do
1478 [ + + ]: 707880 : if (use_gaussian(itypat)) then
1479 : 420048 : sfr=sfr*zion(itypat)
1480 : 420048 : sfi=sfi*zion(itypat)
1481 : : end if
1482 : : else
1483 [ + + ]: 2123640 : fact0=one;if (.not.use_gaussian(itypat)) fact0=one/zion(itypat)
1484 [ + + ]: 6090021 : do ia=ia1,ia2
1485 : : ! Here, take care of the components of m
1486 : 3966381 : fact=fact0*spinat_indx(ispden-1,ia)
1487 : 3966381 : sfr=sfr+phre_ini(ig1,ig2,ig3,ia)*fact
1488 : 6090021 : sfi=sfi-phimag_ini(ig1,ig2,ig3,ia)*fact
1489 : : end do
1490 : : end if
1491 : :
1492 : : ! Charge density integrating to one
1493 [ + + ]: 2831520 : if (use_gaussian(itypat)) then
1494 : 1680192 : rhoat=xnorm*exp(-gsquar*alf2pi2)
1495 : : else
1496 : 1151328 : gmag=sqrt(gsquar)
1497 : 1151328 : jj=1+int(gmag*dqm1)
1498 : 1151328 : diff=gmag-qgrid(jj)
1499 : 1151328 : bb = diff*dqm1
1500 : 1151328 : aa = one-bb
1501 : 1151328 : cc = aa*(aa**2-one)*dq2div6
1502 : 1151328 : dd = bb*(bb**2-one)*dq2div6
1503 [ + + ]: 1151328 : 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 [ + - ]: 536044 : else if (usepaw == 0) then
1507 : : rhoat=(aa*psps%nctab(itypat)%tvalespl(jj,1)+bb*psps%nctab(itypat)%tvalespl(jj+1,1)+&
1508 : 536044 : 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 : 2831520 : rhog(re,ii)=rhog(re,ii)+sfr*rhoat
1516 : 2831520 : 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 : 1120 : ia1=ia2+1
1524 : : end do ! itypat
1525 : :
1526 : : ! Set contribution of unbalanced components to zero
1527 [ + + ]: 496 : 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 : 496 : call fourdp(1,rhog,work,1,mpi_enreg,nfft,1,ngfft,0)
1534 [ + + ]: 6640204 : 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 [ + + + + : 54921 : if (all(abs(spinat(:,:))<tol10)) rhor(:,4)=rhor(:,4)+tol14
+ + + + +
+ ]
1540 : :
1541 : : end if ! nspden==4
1542 : :
1543 [ - + ]: 3096 : ABI_FREE(spinat_indx)
1544 [ + + - + ]: 6192 : ABI_FREE(work)
1545 : :
1546 : : contains
1547 : :
1548 : : !Real and imaginary parts of phase.
1549 : 120879175 : 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 : 120879175 : phr_ini=(x1*x2-y1*y2)*x3-(y1*x2+x1*y2)*y3
1554 : : end function phr_ini
1555 : :
1556 : 120879175 : 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 : 120879175 : phi_ini=(x1*x2-y1*y2)*y3+(y1*x2+x1*y2)*x3
1561 : : end function phi_ini
1562 : :
1563 : 483516700 : function ph1_ini(nri,ig1,ia)
1564 : :
1565 : : real(dp) :: ph1_ini
1566 : : integer,intent(in) :: nri,ig1,ia
1567 : 483516700 : ph1_ini=ph1d(nri,ig1+1+n1+(ia-1)*(2*n1+1))
1568 : 483516700 : end function ph1_ini
1569 : :
1570 : 483516700 : function ph2_ini(nri,ig2,ia)
1571 : :
1572 : : real(dp) :: ph2_ini
1573 : : integer,intent(in) :: nri,ig2,ia
1574 : 483516700 : ph2_ini=ph1d(nri,ig2+1+n2+(ia-1)*(2*n2+1)+natom*(2*n1+1))
1575 : 483516700 : end function ph2_ini
1576 : :
1577 : 483516700 : function ph3_ini(nri,ig3,ia)
1578 : :
1579 : : real(dp) :: ph3_ini
1580 : : integer,intent(in) :: nri,ig3,ia
1581 : 483516700 : ph3_ini=ph1d(nri,ig3+1+n3+(ia-1)*(2*n3+1)+natom*(2*n1+1+2*n2+1))
1582 : 483516700 : end function ph3_ini
1583 : :
1584 : 120879175 : 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 : 120879175 : & ph2_ini(re,ig2,ia),ph2_ini(im,ig2,ia),ph3_ini(re,ig3,ia),ph3_ini(im,ig3,ia))
1590 : 120879175 : end function phre_ini
1591 : :
1592 : 120879175 : 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 : 120879175 : & ph2_ini(re,ig2,ia),ph2_ini(im,ig2,ia),ph3_ini(re,ig3,ia),ph3_ini(im,ig3,ia))
1598 : 120879175 : 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 : 69540 : 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 : 69540 : 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 : 69540 : real(dp),allocatable :: array(:),coord(:,:,:,:),value(:,:,:),integrated(:)
1658 : 69540 : real(dp),allocatable :: value_fft(:,:,:)
1659 : : ! *************************************************************************
1660 : :
1661 [ + + ]: 69540 : if(.not.(present(optrhor))) then
1662 : : ioptrhor=0
1663 : : else
1664 : 16076 : ioptrhor=optrhor
1665 : : end if
1666 : :
1667 [ - + ]: 69540 : 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 [ - + ]: 69540 : 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 : 69540 : nfft_ = nfft
1680 : 69540 : n1=ngfft(1) ; n2=ngfft(2) ; n3=ngfft(3)
1681 : 69540 : 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 [ + + ]: 69540 : if(nspden==2)nitems=5 ! Total density, spin up, spin down, magnetization, zeta
1696 [ + + ]: 69540 : if(nspden==4)nitems=6 ! Total density, x, y, z, magnetization, zeta
1697 : :
1698 [ - + + - : 139080 : ABI_MALLOC(value,(2,2,nitems))
- + - + ]
1699 [ - + + - : 139080 : ABI_MALLOC(iindex,(2,2,nitems))
- + - + ]
1700 [ - + + - : 208620 : ABI_MALLOC(array,(nfft))
+ - + - -
+ - + ]
1701 [ - + + - : 139080 : ABI_MALLOC(integrated,(nitems))
- + - + ]
1702 : :
1703 [ + + ]: 194173 : 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 [ + + ]: 124633 : if(iitems<=nspden)then
1708 [ + + ]: 1858039950 : array(:)=rhor(:,iitems)
1709 : : end if
1710 : : ! Case nspden==2, some computation to be done
1711 [ + + ]: 124633 : if(nspden==2)then
1712 [ + + ]: 39735 : if(iitems==3)then ! Spin down
1713 [ + + ]: 211786279 : array(:)=rhor(:,1)-rhor(:,2)
1714 [ + + ]: 31788 : else if(iitems==4)then ! Magnetization
1715 [ + + ]: 211786279 : array(:)=2*rhor(:,2)-rhor(:,1)
1716 [ + + ]: 23841 : else if(iitems==5)then ! zeta = relative magnetization
1717 : : ! Avoid 0/0: the limit of (x - y) / (x+ y) depends on the direction.
1718 [ + + ]: 211786279 : array(:)=zero
1719 [ + + + + ]: 211786279 : 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 [ + + ]: 84898 : else if(nspden==4)then
1723 [ + + ]: 27966 : if(iitems==5)then ! Magnetization
1724 [ + + ]: 46792561 : array(:)=sqrt(rhor(:,2)**2+rhor(:,3)**2+rhor(:,4)**2)
1725 [ + + ]: 23305 : else if(iitems==6)then ! zeta = relative magnetization
1726 [ + + ]: 46792561 : 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 [ + + ]: 2586983909 : do ifft = 1, nfft_
1732 [ + + ]: 2586983909 : 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 [ + + ]: 373899 : do indsign=1,2 ! Find alternatively the maximum and the minimum
1740 : 249266 : iisign=3-2*indsign
1741 : :
1742 [ + - ]: 249266 : if (nfft_ > 1) then
1743 : : ! Initialize the two first values
1744 : 249266 : value1=array(istart + 1) ; value2=array(istart + 2)
1745 : 249266 : index1=1 ; index2=2
1746 : :
1747 : : ! Ordering, if needed
1748 [ + + ]: 249266 : if( iisign*(value2+tol12) > iisign*(value1)) then
1749 : 124633 : temp=value2 ; value2=value1 ; value1=temp
1750 : 124633 : index1=2 ; index2=1
1751 : : end if
1752 : :
1753 : : ! Integration, if relevant
1754 [ + + ]: 249266 : if(present(ucvol).and. indsign==1)then
1755 : 124633 : 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 [ + + ]: 5173469286 : do ifft = 3, nfft_
1768 : :
1769 : 5173220020 : temp=array(istart + ifft)
1770 [ + + ]: 5173220020 : if(present(ucvol).and. indsign==1)integrated(iitems) = integrated(iitems)+temp
1771 : : ! Compares it to the second value
1772 [ + + ]: 5173469286 : if( iisign*(temp+tol12) > iisign*value2 ) then
1773 : : ! Compare it to the first value
1774 [ + + ]: 215798913 : if( iisign*(temp+tol12) > iisign*value1 ) then
1775 : 213320312 : value2=value1 ; index2=index1
1776 : 213320312 : value1=temp ; index1=ifft
1777 : : else
1778 : : value2=temp ; index2=ifft
1779 : : end if
1780 : : end if
1781 : :
1782 : : end do ! ifft
1783 : :
1784 : 249266 : value(1,indsign,iitems)=value1
1785 : 249266 : value(2,indsign,iitems)=value2
1786 : 249266 : iindex(1,indsign,iitems)=index1
1787 : 373899 : 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 [ + - ]: 194173 : if(present(ucvol))then
1797 : 124633 : nfftot=ngfft(1) * ngfft(2) * ngfft(3)
1798 : 124633 : 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 [ - + ]: 69540 : 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 : 69540 : spaceComm=mpi_enreg%comm_fft; reduce=.false.
1820 [ + + ]: 69540 : if(mpi_enreg%nproc_fft>1) then
1821 : 4304 : spaceComm=mpi_enreg%comm_fft; reduce=.true.
1822 [ - + ]: 65236 : else if(mpi_enreg%nproc_wvl>1) then
1823 : 0 : spaceComm=mpi_enreg%comm_wvl; reduce=.true.
1824 : : end if
1825 : 69540 : nproc=xmpi_comm_size(spaceComm)
1826 : 69540 : me=xmpi_comm_rank(spaceComm)
1827 : :
1828 [ + + ]: 69540 : if (reduce) then
1829 : :
1830 : : ! Communicate all data to all processors with only two global communications
1831 [ - + + - : 17216 : ABI_MALLOC(value_fft,(5,nitems,nproc))
- + + - +
- + - - +
- + ]
1832 [ - + + - : 17216 : ABI_MALLOC(index_fft,(2,2,nitems,nproc))
- + + - +
- + - - +
- + ]
1833 [ + + + + : 74336 : value_fft(:,:,:)=zero
+ + ]
1834 [ + + + + : 84504 : index_fft(:,:,:,:)=0
+ + + + ]
1835 [ + + ]: 9180 : value_fft(1,:,me + 1)=value(1,1,:)
1836 [ + + ]: 9180 : value_fft(2,:,me + 1)=value(2,1,:)
1837 [ + + ]: 9180 : value_fft(3,:,me + 1)=value(1,2,:)
1838 [ + + ]: 9180 : value_fft(4,:,me + 1)=value(2,2,:)
1839 [ + - + + ]: 9180 : if(present(ucvol))value_fft(5,:,me + 1)=integrated(:)
1840 [ + + + + : 38436 : index_fft(:,:,:,me + 1)=iindex(:,:,:)
+ + ]
1841 : 4304 : call xmpi_sum(value_fft,spaceComm,ierr)
1842 : 4304 : 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 [ + + ]: 9180 : do iitems=1,nitems
1847 : :
1848 [ + - + + ]: 15044 : if(present(ucvol))integrated(iitems)=sum(value_fft(5,iitems,1:nproc))
1849 : :
1850 [ + + ]: 18932 : do indsign=1,2 ! Find alternatively the maximum and the minimum
1851 : 9752 : iisign=3-2*indsign
1852 : :
1853 : : ! Initialisation
1854 : 9752 : value1=value_fft(2*indsign-1,iitems,1)
1855 : 9752 : value2=value_fft(2*indsign ,iitems,1)
1856 : 9752 : index1=index_fft(1,indsign,iitems,1)
1857 : 9752 : index2=index_fft(2,indsign,iitems,1)
1858 : :
1859 : : ! Loop
1860 [ + + ]: 30088 : do iproc=1, nproc, 1
1861 [ + + ]: 70760 : do ii=1,2
1862 [ + + ]: 61008 : if(iproc>1 .or. ii==2)then
1863 : :
1864 : 30920 : temp=value_fft(ii+2*(indsign-1),iitems,iproc)
1865 : : ! Compares it to the second value
1866 [ + + ]: 30920 : if( iisign*(temp+tol12) > iisign*value2 ) then
1867 : : ! Compare it to the first value
1868 [ + + ]: 13340 : if( iisign*(temp+tol12) > iisign*value1 ) then
1869 : 6488 : value2=value1 ; index2=index1
1870 : 6488 : value1=temp ; index1=index_fft(ii,indsign,iitems,iproc)
1871 : : else
1872 : 6852 : 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 : 9752 : value(1,indsign,iitems)=value1
1881 : 9752 : value(2,indsign,iitems)=value2
1882 : 9752 : iindex(1,indsign,iitems)=index1
1883 : 14628 : iindex(2,indsign,iitems)=index2
1884 : :
1885 : : end do ! iisign
1886 : : end do ! iitems
1887 : :
1888 [ - + ]: 4304 : ABI_FREE(value_fft)
1889 [ - + ]: 8608 : 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 [ - + + - : 139080 : ABI_MALLOC(coord,(3,2,2,nitems))
- + - + ]
1897 [ + + ]: 194173 : do iitems=1,nitems
1898 [ + + ]: 443439 : do indsign=1,2
1899 [ + + ]: 872431 : do ii=1,2
1900 : 498532 : index1=iindex(ii,indsign,iitems)
1901 : 498532 : i3=(index1-1)/n1/n2
1902 : 498532 : i2=(index1-1-i3*n1*n2)/n1
1903 : 498532 : i1=index1-1-i3*n1*n2-i2*n1
1904 : 498532 : coord(1,ii,indsign,iitems)=dble(i1)/dble(n1)+tol12
1905 : 498532 : coord(2,ii,indsign,iitems)=dble(i2)/dble(n2)+tol12
1906 : 747798 : 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 [ + + + + ]: 69540 : if (mpi_enreg%paral_kgb==0.or.mpi_enreg%me_fft==0) then
1918 : : if(.true.)then
1919 [ + + ]: 189473 : do iitems=1,nitems
1920 : :
1921 [ + + ]: 122139 : if(ioptrhor==4 .and. iitems>2)exit
1922 : :
1923 : 120527 : select case (ioptrhor)
1924 : : case(0)
1925 : :
1926 [ + + ]: 120527 : if(iitems==1) write(message,'(a)')' Total charge density [el/Bohr^3]'
1927 [ + + ]: 120527 : if(nspden==2)then
1928 [ + + ]: 39000 : if(iitems==2) write(message,'(a)')' Spin up density [el/Bohr^3]'
1929 [ + + ]: 39000 : if(iitems==3) write(message,'(a)')' Spin down density [el/Bohr^3]'
1930 [ + + ]: 39000 : if(iitems==4) write(message,'(a)')' Magnetization (spin up - spin down) [el/Bohr^3]'
1931 [ + + ]: 39000 : if(iitems==5) write(message,'(a)')' Relative magnetization (=zeta, between -1 and 1) '
1932 [ + + ]: 81527 : else if(nspden==4)then
1933 [ + + ]: 27450 : if(iitems==2) write(message,'(a)')' x component of magnetization [el/Bohr^3]'
1934 [ + + ]: 27450 : if(iitems==3) write(message,'(a)')' y component of magnetization [el/Bohr^3]'
1935 [ + + ]: 27450 : if(iitems==4) write(message,'(a)')' z component of magnetization [el/Bohr^3]'
1936 [ + + ]: 27450 : if(iitems==5) write(message,'(a)')' Magnetization (absolute value) [el/Bohr^3]'
1937 [ + + ]: 27450 : 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 [ + + + + : 122143 : 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 : 122137 : call wrtout(iout,message,'COLL')
2006 : :
2007 : 122137 : write(message,'(a,es13.4,a,3f10.4)') ') Maximum= ',&
2008 : 244274 : & value(1,1,iitems),' at reduced coord.',coord(:,1,1,iitems)
2009 : 122137 : call wrtout(iout,message,'COLL')
2010 [ + + ]: 122137 : if(option==2)then
2011 : 9432 : write(message,'(a,es13.4,a,3f10.4)') ')Next maximum= ',&
2012 : 18864 : & value(2,1,iitems),' at reduced coord.',coord(:,2,1,iitems)
2013 : 9432 : call wrtout(iout,message,'COLL')
2014 : : end if
2015 : 122137 : write(message,'(a,es13.4,a,3f10.4)') ') Minimum= ',&
2016 : 244274 : & value(1,2,iitems),' at reduced coord.',coord(:,1,2,iitems)
2017 : 122137 : call wrtout(iout,message,'COLL')
2018 [ + + ]: 122137 : if(option==2)then
2019 : 9432 : write(message,'(a,es13.4,a,3f10.4)') ')Next minimum= ',&
2020 : 18864 : & value(2,2,iitems),' at reduced coord.',coord(:,2,2,iitems)
2021 : 9432 : call wrtout(iout,message,'COLL')
2022 : : end if
2023 [ + - ]: 189473 : if(present(ucvol))then
2024 [ + + + + : 122137 : if(.not.(nspden==2.and.iitems==5) .and. .not.(nspden==4.and.iitems==6))then
+ + + + ]
2025 [ + + ]: 109595 : if(abs(integrated(iitems))<tol10)integrated(iitems)=zero
2026 : 109595 : write(message,'(a,es13.4)')' Integrated= ',integrated(iitems)
2027 : 109595 : 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 [ - + ]: 69540 : ABI_FREE(coord)
2227 [ - + ]: 69540 : ABI_FREE(value)
2228 [ - + ]: 69540 : ABI_FREE(iindex)
2229 [ - + ]: 69540 : ABI_FREE(integrated)
2230 : :
2231 [ - + - + ]: 69540 : 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 : : !!***
|
