Line data Source code
1 : !!****m* ABINIT/m_pawrhoij
2 : !! NAME
3 : !! m_pawrhoij
4 : !!
5 : !! FUNCTION
6 : !! This module contains the definition of the pawrhoij_type structured datatype,
7 : !! as well as related functions and methods.
8 : !! pawrhoij_type variables define rhoij occupancies matrixes used within PAW formalism.
9 : !!
10 : !! COPYRIGHT
11 : !! Copyright (C) 2012-2026 ABINIT group (MT, FJ)
12 : !! This file is distributed under the terms of the
13 : !! GNU General Public License, see ~abinit/COPYING
14 : !! or http://www.gnu.org/copyleft/gpl.txt .
15 : !!
16 : !! NOTES
17 : !! FOR DEVELOPPERS: in order to preserve the portability of libPAW library,
18 : !! please consult ~abinit/src/??_libpaw/libpaw-coding-rules.txt
19 : !!
20 : !! SOURCE
21 :
22 : #include "libpaw.h"
23 :
24 : MODULE m_pawrhoij
25 :
26 : USE_DEFS
27 : USE_MSG_HANDLING
28 : USE_MPI_WRAPPERS
29 : USE_MEMORY_PROFILING
30 : #ifdef LIBPAW_HAVE_NETCDF
31 : use netcdf
32 : #endif
33 :
34 : use m_libpaw_tools, only : libpaw_flush, libpaw_to_upper
35 :
36 : use m_paw_io, only : pawio_print_ij
37 : use m_pawang, only : pawang_type
38 : use m_pawtab, only : pawtab_type
39 : use m_paral_atom, only : get_my_atmtab, free_my_atmtab, get_my_natom
40 :
41 : implicit none
42 :
43 : private
44 :
45 : !public procedures.
46 : public :: pawrhoij_alloc
47 : public :: pawrhoij_free
48 : public :: pawrhoij_nullify
49 : public :: pawrhoij_copy
50 : public :: pawrhoij_gather
51 : public :: pawrhoij_bcast
52 : public :: pawrhoij_redistribute
53 : public :: pawrhoij_io
54 : public :: pawrhoij_unpack
55 : public :: pawrhoij_init_unpacked
56 : public :: pawrhoij_free_unpacked
57 : public :: pawrhoij_filter
58 : public :: pawrhoij_inquire_dim
59 : public :: pawrhoij_print_rhoij
60 : public :: pawrhoij_symrhoij
61 :
62 : public :: pawrhoij_mpisum_unpacked
63 : interface pawrhoij_mpisum_unpacked
64 : module procedure pawrhoij_mpisum_unpacked_1D
65 : module procedure pawrhoij_mpisum_unpacked_2D
66 : end interface pawrhoij_mpisum_unpacked
67 :
68 : !private procedures.
69 : private :: pawrhoij_isendreceive_getbuffer
70 : private :: pawrhoij_isendreceive_fillbuffer
71 : !!***
72 :
73 : !!****t* m_pawrhoij/pawrhoij_type
74 : !! NAME
75 : !! pawrhoij_type
76 : !!
77 : !! FUNCTION
78 : !! This structured datatype contains rhoij quantities (occucpancies)
79 : !! and related data, used in PAW calculations.
80 : !!
81 : !! SOURCE
82 :
83 : type,public :: pawrhoij_type
84 :
85 : !Integer scalars
86 :
87 : integer :: cplex_rhoij
88 : ! cplex_rhoij=1 if rhoij are real
89 : ! cplex_rhoij=2 if rhoij are complex (spin-orbit, pawcpxocc=2, ...)
90 :
91 : integer :: itypat
92 : ! itypat=type of the atom
93 :
94 : integer :: lmn_size
95 : ! Number of (l,m,n) elements for the paw basis
96 :
97 : integer :: lmn2_size
98 : ! lmn2_size=lmn_size*(lmn_size+1)/2
99 : ! where lmn_size is the number of (l,m,n) elements for the paw basis
100 :
101 : integer :: lmnmix_sz=0
102 : ! lmnmix_sz=number of (lmn,lmn_prime) verifying l<=lmix and l_prime<=lmix
103 : ! i.e. number of rhoij elements being mixed during SCF cycle
104 : ! lmnmix_sz=0 if mixing data are not used
105 :
106 : integer :: ngrhoij=0
107 : ! First dimension of array grhoij
108 :
109 : integer :: nrhoijsel=0
110 : ! nrhoijsel
111 : ! Number of non-zero values of rhoij
112 : ! This is the size of rhoijp(:,:) (see below in this datastructure)
113 :
114 : integer :: nspden
115 : ! Number of spin-density components for rhoij (may be different from nspden for density)
116 :
117 : integer :: nspinor
118 : ! Number of spinorial components
119 :
120 : integer :: nsppol
121 : ! Number of independent spin-components
122 :
123 : integer :: qphase
124 : ! qphase=2 if rhoij contain a exp(-i.q.r) phase (as in the q<>0 RF case), 1 if not
125 : ! (this may change the ij symmetry)
126 :
127 : integer :: use_rhoij_=0
128 : ! 1 if pawrhoij%rhoij_ is allocated
129 :
130 : integer :: use_rhoijp=0
131 : ! 1 if pawrhoij%rhoijp and pawrhoij%rhoijselect are allocated
132 :
133 : integer :: use_rhoijres=0
134 : ! 1 if pawrhoij%rhoijres is allocated
135 :
136 : !Integer arrays
137 :
138 : integer, allocatable :: kpawmix(:)
139 : ! kpawmix(lmnmix_sz)
140 : ! Indirect array selecting the elements of rhoij
141 : ! being mixed during SCF cycle
142 :
143 : integer, allocatable :: rhoijselect(:)
144 : ! rhoijselect(lmn2_size)
145 : ! Indirect array selecting the non-zero elements of rhoij:
146 : ! rhoijselect(isel,ispden)=klmn if rhoij(klmn,ispden) is non-zero
147 :
148 : !Real (real(dp)) arrays
149 :
150 : real(dp), allocatable :: grhoij (:,:,:)
151 : ! grhoij(ngrhoij,cplex_rhoij*qphase*lmn2_size,nspden)
152 : ! Gradients of Rho_ij wrt xred, strains, ... (non-packed storage)
153 :
154 : real(dp), allocatable :: rhoij_ (:,:)
155 : ! rhoij_(cplex_rhoij*qphase*lmn2_size,nspden)
156 : ! Array used to (temporary) store Rho_ij in a non-packed storage mode
157 :
158 : real(dp), allocatable :: rhoijp (:,:)
159 : ! rhoijp(cplex_rhoij*qphase*lmn2_size,nspden)
160 : ! Augmentation waves occupancies Rho_ij in PACKED STORAGE (only non-zero elements are stored)
161 :
162 : real(dp), allocatable :: rhoijres (:,:)
163 : ! rhoijres(cplex_rhoij*qphase*lmn2_size,nspden)
164 : ! Rho_ij residuals during SCF cycle (non-packed storage)
165 :
166 : ! ==== Storage for the 1st dimension ====
167 : ! For each klmn=ij:
168 : ! When RHOij is complex (cplex=2):
169 : ! rhoij(2*ij-1,:) contains the real part
170 : ! rhoij(2*ij ,:) contains the imaginary part
171 : ! When a exp(-i.q.r) phase is included (qphase=2):
172 : ! rhoij(1:cplex_dij*lmn2_size,:)
173 : ! contains the real part of the phase, i.e. RHO_ij*cos(q.r)
174 : ! rhoij(cplex_dij*lmn2_size+1:2*cplex_dij*lmn2_size,:)
175 : ! contains the imaginary part of the phase, i.e. RHO_ij*sin(q.r)
176 : ! ==== Storage for the 2nd dimension ====
177 : ! No magnetism
178 : ! rhoij(:,1) contains rhoij
179 : ! Collinear magnetism
180 : ! rhoij(:,1) contains rhoij^up
181 : ! rhoij(:,2) contains rhoij^dowm
182 : ! Non-collinear magnetism
183 : ! rhoij(:,1) contains rhoij
184 : ! rhoij(:,2) contains rhoij magnetization along x
185 : ! rhoij(:,3) contains rhoij magnetization along y
186 : ! rhoij(:,4) contains rhoij magnetization along z
187 :
188 : end type pawrhoij_type
189 : !!***
190 :
191 : CONTAINS
192 :
193 : !===========================================================
194 : !!***
195 :
196 : !----------------------------------------------------------------------
197 :
198 : !!****f* m_pawrhoij/pawrhoij_alloc
199 : !! NAME
200 : !! pawrhoij_alloc
201 : !!
202 : !! FUNCTION
203 : !! Initialize and allocate a pawrhoij datastructure
204 : !!
205 : !! INPUTS
206 : !! [comm_atom] = communicator over atoms (OPTIONAL)
207 : !! cplex_rhoij=1 if rhoij are real, 2 if rhoij are complex (spin-orbit, pawcpxocc=2, ...)
208 : !! [my_atmtab(:)] = Index of atoms treated by current proc (OPTIONAL)
209 : !! nspden=number of spin-components for rhoij
210 : !! nsppol=number of spinorial components for rhoij
211 : !! nsppol=number of independant spin-components for rhoij
212 : !! typat(:)=types of atoms
213 : !! [lmnsize(:)]=array of (l,m,n) sizes for rhoij for each type of atom (OPTIONAL)
214 : !! must be present if [pawtab] argument is not passed
215 : !! [ngrhoij]=number of gradients to be allocated (OPTIONAL, default=0)
216 : !! [nlmnmix]=number of rhoij elements to be mixed during SCF cycle (OPTIONAL, default=0)
217 : !! [pawtab(:)] <type(pawtab_type)>=paw tabulated starting data (OPTIONAL)
218 : !! must be present if [lmnsize(:)] argument is not passed
219 : !! [qphase]=2 if the rhoij contain a exp(iqR) phase, 1 otherwise (OPTIONAL, default=1)
220 : !! Typical use: 1st-order rhoij at q<>0
221 : !! [use_rhoij_]=1 if pawrhoij(:)%rhoij_ has to be allocated (OPTIONAL, default=0)
222 : !! [use_rhoijp]=1 if pawrhoij(:)%rhoijp has to be allocated (OPTIONAL, default=1)
223 : !! (in that case, pawrhoij%rhoijselect is also allocated)
224 : !! [use_rhoijres]=1 if pawrhoij(:)%rhoijres has to be allocated (OPTIONAL, default=0)
225 : !!
226 : !! SIDE EFFECTS
227 : !! pawrhoij(:)<type(pawrhoij_type)>= rhoij datastructure
228 : !!
229 : !! NOTES
230 : !! One of the two optional arguments lmnsize(:) or pawtab(:) must be present !
231 : !! If both are present, only pawtab(:) is used.
232 : !!
233 : !! SOURCE
234 :
235 19996 : subroutine pawrhoij_alloc(pawrhoij,cplex_rhoij,nspden,nspinor,nsppol,typat,&
236 19996 : & lmnsize,ngrhoij,nlmnmix,pawtab,qphase,use_rhoij_,use_rhoijp,& ! Optional
237 19996 : & use_rhoijres,comm_atom,mpi_atmtab) ! Optional
238 :
239 : !Arguments ------------------------------------
240 : !scalars
241 : integer,intent(in) :: cplex_rhoij,nspden,nspinor,nsppol
242 : integer,optional,intent(in):: comm_atom,ngrhoij,nlmnmix,qphase
243 : integer,optional,intent(in):: use_rhoij_,use_rhoijp,use_rhoijres
244 : integer,optional,target,intent(in) :: mpi_atmtab(:)
245 : !arrays
246 : integer,intent(in) :: typat(:)
247 : integer,optional,target,intent(in) :: lmnsize(:)
248 : type(pawrhoij_type),intent(inout) :: pawrhoij(:)
249 : type(pawtab_type),optional,intent(in) :: pawtab(:)
250 :
251 : !Local variables-------------------------------
252 : !scalars
253 : integer :: irhoij,irhoij_,itypat,lmn2_size,my_comm_atom,my_qphase, nn1,natom,nrhoij
254 : logical :: has_rhoijp,my_atmtab_allocated,paral_atom
255 : character(len=500) :: msg
256 : !array
257 19996 : integer,pointer :: lmn_size(:),my_atmtab(:)
258 :
259 : ! *************************************************************************
260 :
261 19996 : nrhoij=size(pawrhoij);natom=size(typat)
262 19996 : if (nrhoij>natom) then
263 0 : msg=' wrong sizes (1) !'
264 0 : LIBPAW_BUG(msg)
265 : end if
266 :
267 : !Select lmn_size for each atom type
268 19996 : if (present(pawtab)) then
269 15009 : nn1=size(pawtab)
270 64221 : if (maxval(typat)>nn1) then
271 0 : msg=' wrong sizes (2) !'
272 0 : LIBPAW_BUG(msg)
273 : end if
274 45027 : LIBPAW_POINTER_ALLOCATE(lmn_size,(nn1))
275 41057 : do itypat=1,nn1
276 41057 : lmn_size(itypat)=pawtab(itypat)%lmn_size
277 : end do
278 4987 : else if (present(lmnsize)) then
279 4987 : nn1=size(lmnsize)
280 15483 : if (maxval(typat)>nn1) then
281 0 : msg=' wrong sizes (3) !'
282 0 : LIBPAW_BUG(msg)
283 : end if
284 4987 : lmn_size => lmnsize
285 : else
286 0 : msg=' one of the 2 arguments pawtab or lmnsize must be present !'
287 0 : LIBPAW_BUG(msg)
288 : end if
289 :
290 : !Set up parallelism over atoms
291 19996 : paral_atom=(present(comm_atom).and.(nrhoij/=natom))
292 19996 : nullify(my_atmtab);if (present(mpi_atmtab)) my_atmtab => mpi_atmtab
293 19996 : my_comm_atom=xmpi_comm_self;if (present(comm_atom)) my_comm_atom=comm_atom
294 19996 : call get_my_atmtab(my_comm_atom,my_atmtab,my_atmtab_allocated,paral_atom,natom)
295 :
296 19996 : my_qphase=1;if (present(qphase)) my_qphase=qphase
297 :
298 19996 : if (nrhoij>0) then
299 72913 : do irhoij=1,nrhoij
300 52951 : irhoij_=irhoij;if (paral_atom) irhoij_=my_atmtab(irhoij)
301 52951 : itypat=typat(irhoij_)
302 :
303 52951 : lmn2_size=lmn_size(itypat)*(lmn_size(itypat)+1)/2
304 :
305 : ! Scalars initializations
306 52951 : pawrhoij(irhoij)%cplex_rhoij=cplex_rhoij
307 52951 : pawrhoij(irhoij)%qphase=my_qphase
308 52951 : pawrhoij(irhoij)%itypat=itypat
309 52951 : pawrhoij(irhoij)%lmn_size=lmn_size(itypat)
310 52951 : pawrhoij(irhoij)%lmn2_size=lmn2_size
311 52951 : pawrhoij(irhoij)%nspden=nspden
312 52951 : pawrhoij(irhoij)%nspinor=nspinor
313 52951 : pawrhoij(irhoij)%nsppol=nsppol
314 52951 : pawrhoij(irhoij)%nrhoijsel=0
315 52951 : pawrhoij(irhoij)%lmnmix_sz=0
316 52951 : pawrhoij(irhoij)%ngrhoij=0
317 52951 : pawrhoij(irhoij)%use_rhoij_=0
318 52951 : pawrhoij(irhoij)%use_rhoijres=0
319 :
320 : ! Arrays allocations
321 52951 : has_rhoijp=.true.; if (present(use_rhoijp)) has_rhoijp=(use_rhoijp>0)
322 33288 : if (has_rhoijp) then
323 27367 : pawrhoij(irhoij)%use_rhoijp=1
324 82101 : LIBPAW_ALLOCATE(pawrhoij(irhoij)%rhoijselect,(lmn2_size))
325 109468 : LIBPAW_ALLOCATE(pawrhoij(irhoij)%rhoijp,(cplex_rhoij*my_qphase*lmn2_size,nspden))
326 1386803 : pawrhoij(irhoij)%rhoijselect(:)=0
327 2308163 : pawrhoij(irhoij)%rhoijp(:,:)=zero
328 : end if
329 :
330 52951 : if (present(ngrhoij)) then
331 4856 : if (ngrhoij>0) then
332 0 : pawrhoij(irhoij)%ngrhoij=ngrhoij
333 0 : LIBPAW_ALLOCATE(pawrhoij(irhoij)%grhoij,(ngrhoij,cplex_rhoij*my_qphase*lmn2_size,nspden))
334 0 : pawrhoij(irhoij)%grhoij=zero
335 : end if
336 : end if
337 52951 : if (present(nlmnmix)) then
338 4856 : if (nlmnmix>0) then
339 1374 : pawrhoij(irhoij)%lmnmix_sz=nlmnmix
340 4122 : LIBPAW_ALLOCATE(pawrhoij(irhoij)%kpawmix,(nlmnmix))
341 71077 : pawrhoij(irhoij)%kpawmix=0
342 : end if
343 : end if
344 52951 : if (present(use_rhoij_)) then
345 24411 : if (use_rhoij_>0) then
346 13516 : pawrhoij(irhoij)%use_rhoij_=use_rhoij_
347 54064 : LIBPAW_ALLOCATE(pawrhoij(irhoij)%rhoij_,(cplex_rhoij*my_qphase*lmn2_size,nspden))
348 532556 : pawrhoij(irhoij)%rhoij_=zero
349 : end if
350 : end if
351 72913 : if (present(use_rhoijres)) then
352 4856 : if (use_rhoijres>0) then
353 1374 : pawrhoij(irhoij)%use_rhoijres=use_rhoijres
354 5496 : LIBPAW_ALLOCATE(pawrhoij(irhoij)%rhoijres,(cplex_rhoij*my_qphase*lmn2_size,nspden))
355 91904 : pawrhoij(irhoij)%rhoijres=zero
356 : end if
357 : end if
358 :
359 : end do
360 : end if
361 :
362 19996 : if (present(pawtab)) then
363 15009 : LIBPAW_POINTER_DEALLOCATE(lmn_size)
364 : end if
365 :
366 : !Destroy atom table used for parallelism
367 19996 : call free_my_atmtab(my_atmtab,my_atmtab_allocated)
368 :
369 19996 : end subroutine pawrhoij_alloc
370 : !!***
371 :
372 : !----------------------------------------------------------------------
373 :
374 : !!****f* m_pawrhoij/pawrhoij_free
375 : !! NAME
376 : !! pawrhoij_free
377 : !!
378 : !! FUNCTION
379 : !! Destroy a pawrhoij datastructure
380 : !!
381 : !! SIDE EFFECTS
382 : !! pawrhoij(:)<type(pawrhoij_type)>= rhoij datastructure
383 : !!
384 : !! SOURCE
385 :
386 25346 : subroutine pawrhoij_free(pawrhoij)
387 :
388 : !Arguments ------------------------------------
389 : !arrays
390 : type(pawrhoij_type),intent(inout) :: pawrhoij(:)
391 :
392 : !Local variables-------------------------------
393 : !scalars
394 : integer :: irhoij,nrhoij
395 :
396 : ! *************************************************************************
397 :
398 25346 : nrhoij=size(pawrhoij)
399 :
400 25346 : if (nrhoij>0) then
401 88085 : do irhoij=1,nrhoij
402 64307 : pawrhoij(irhoij)%cplex_rhoij=1
403 64307 : pawrhoij(irhoij)%qphase=1
404 64307 : pawrhoij(irhoij)%nrhoijsel=0
405 64307 : pawrhoij(irhoij)%ngrhoij=0
406 64307 : pawrhoij(irhoij)%lmnmix_sz=0
407 64307 : pawrhoij(irhoij)%use_rhoij_=0
408 64307 : pawrhoij(irhoij)%use_rhoijp=0
409 64307 : pawrhoij(irhoij)%use_rhoijres=0
410 64307 : if (allocated(pawrhoij(irhoij)%rhoijp)) then
411 33063 : LIBPAW_DEALLOCATE(pawrhoij(irhoij)%rhoijp)
412 : end if
413 64307 : if (allocated(pawrhoij(irhoij)%rhoijselect)) then
414 33063 : LIBPAW_DEALLOCATE(pawrhoij(irhoij)%rhoijselect)
415 : end if
416 64307 : if (allocated(pawrhoij(irhoij)%grhoij)) then
417 0 : LIBPAW_DEALLOCATE(pawrhoij(irhoij)%grhoij)
418 : end if
419 64307 : if (allocated(pawrhoij(irhoij)%kpawmix)) then
420 2944 : LIBPAW_DEALLOCATE(pawrhoij(irhoij)%kpawmix)
421 : end if
422 64307 : if (allocated(pawrhoij(irhoij)%rhoij_)) then
423 8 : LIBPAW_DEALLOCATE(pawrhoij(irhoij)%rhoij_)
424 : end if
425 88085 : if (allocated(pawrhoij(irhoij)%rhoijres)) then
426 2944 : LIBPAW_DEALLOCATE(pawrhoij(irhoij)%rhoijres)
427 : end if
428 : end do
429 : end if
430 :
431 25346 : end subroutine pawrhoij_free
432 : !!***
433 :
434 : !----------------------------------------------------------------------
435 :
436 : !!****f* m_pawrhoij/pawrhoij_nullify
437 : !! NAME
438 : !! pawrhoij_nullify
439 : !!
440 : !! FUNCTION
441 : !! Nullify (initialize to null) a pawrhoij datastructure
442 : !!
443 : !! SIDE EFFECTS
444 : !! pawrhoij(:)<type(pawrhoij_type)>= rhoij datastructure
445 : !!
446 : !! SOURCE
447 :
448 3595 : subroutine pawrhoij_nullify(pawrhoij)
449 :
450 : !Arguments ------------------------------------
451 : !arrays
452 : type(pawrhoij_type),intent(inout) :: pawrhoij(:)
453 :
454 : !Local variables-------------------------------
455 : !scalars
456 : integer :: irhoij,nrhoij
457 :
458 : ! *************************************************************************
459 :
460 : ! MGPAW: This one could be removed/renamed,
461 : ! variables can be initialized in the datatype declaration
462 : ! Do we need to expose this in the public API?
463 :
464 3595 : nrhoij=size(pawrhoij)
465 :
466 3595 : if (nrhoij>0) then
467 7603 : do irhoij=1,nrhoij
468 4773 : pawrhoij(irhoij)%cplex_rhoij=1
469 4773 : pawrhoij(irhoij)%qphase=1
470 4773 : pawrhoij(irhoij)%nrhoijsel=0
471 4773 : pawrhoij(irhoij)%ngrhoij=0
472 4773 : pawrhoij(irhoij)%lmnmix_sz=0
473 4773 : pawrhoij(irhoij)%use_rhoij_=0
474 4773 : pawrhoij(irhoij)%use_rhoijp=0
475 7603 : pawrhoij(irhoij)%use_rhoijres=0
476 : end do
477 : end if
478 :
479 3595 : end subroutine pawrhoij_nullify
480 : !!***
481 :
482 : !----------------------------------------------------------------------
483 :
484 : !!****f* m_pawrhoij/pawrhoij_copy
485 : !! NAME
486 : !! pawrhoij_copy
487 : !!
488 : !! FUNCTION
489 : !! Copy one pawrhoij datastructure into another
490 : !! Can take into accound changes of dimensions
491 : !! Can copy a shared pawrhoij into distributed ones (when parallelism is activated)
492 : !!
493 : !! INPUTS
494 : !! keep_cplex= optional argument (logical, default=.TRUE.)
495 : !! if .TRUE. pawrhoij_out(:)%cplex_rhoij is NOT MODIFIED,
496 : !! even if different from pawrhoij_in(:)%cplex_rhoij
497 : !! keep_qphase= optional argument (logical, default=.TRUE.)
498 : !! if .TRUE. pawrhoij_out(:)%cplex_rhoij is NOT MODIFIED,
499 : !! even if different from pawrhoij_in(:)%qphase
500 : !! keep_itypat= optional argument (logical, default=.FALSE.)
501 : !! if .TRUE. pawrhoij_out(:)%ityp is NOT MODIFIED,
502 : !! even if different from pawrhoij_in(:)%ityp
503 : !! keep_nspden= optional argument (logical, default=.TRUE.)
504 : !! if .TRUE. pawrhoij_out(:)%nspden is NOT MODIFIED,
505 : !! even if different from pawrhoij_in(:)%nspden
506 : !! mpi_atmtab(:)=--optional-- indexes of the atoms treated by current proc
507 : !! comm_atom=--optional-- MPI communicator over atoms
508 : !! pawrhoij_in(:)<type(pawrhoij_type)>= input rhoij datastructure
509 : !!
510 : !! SIDE EFFECTS
511 : !! pawrhoij_out(:)<type(pawrhoij_type)>= output rhoij datastructure
512 : !!
513 : !! NOTES
514 : !! In case of a single copy operation pawrhoij_out must have been allocated.
515 : !!
516 : !! SOURCE
517 :
518 8947 : subroutine pawrhoij_copy(pawrhoij_in,pawrhoij_cpy, &
519 : & keep_cplex,keep_qphase,keep_itypat,keep_nspden,& ! optional arguments
520 8947 : & mpi_atmtab,comm_atom) ! optional arguments (parallelism)
521 :
522 : !Arguments ------------------------------------
523 : !scalars
524 : integer,optional,intent(in) :: comm_atom
525 : logical,intent(in),optional :: keep_cplex,keep_qphase,keep_itypat,keep_nspden
526 : !arrays
527 : integer,optional,target,intent(in) :: mpi_atmtab(:)
528 : type(pawrhoij_type),intent(in) :: pawrhoij_in(:)
529 : type(pawrhoij_type),intent(inout),target :: pawrhoij_cpy(:)
530 :
531 : !Local variables-------------------------------
532 : !scalars
533 : integer :: cplex,cplex_in,cplex_out,i_in,i_out,ilmn,iphase
534 : integer :: irhoij,ispden,jrhoij,lmn2_size_in,lmn2_size_out,lmnmix,my_comm_atom,my_nrhoij
535 : integer :: ngrhoij,nrhoij_in,nrhoij_max,nrhoij_out,nselect,nselect_out
536 : integer :: nspden_in,nspden_out,paral_case,qphase,qphase_in,qphase_out
537 : integer :: use_rhoij_,use_rhoijp,use_rhoijres
538 : logical :: change_dim,keep_cplex_,keep_qphase_,keep_itypat_,keep_nspden_,my_atmtab_allocated,paral_atom
539 : character(len=500) :: msg
540 : !arrays
541 8947 : integer,pointer :: my_atmtab(:)
542 8947 : integer,allocatable :: nlmn(:),typat(:)
543 8947 : type(pawrhoij_type),pointer :: pawrhoij_out(:)
544 :
545 : ! *************************************************************************
546 :
547 : !Retrieve sizes
548 8947 : nrhoij_in=size(pawrhoij_in);nrhoij_out=size(pawrhoij_cpy)
549 :
550 : !Init flags
551 8947 : keep_cplex_=.true.
552 60 : if (present(keep_cplex)) keep_cplex_=keep_cplex
553 8947 : keep_qphase_=.true.
554 8947 : if (present(keep_qphase)) keep_qphase_=keep_qphase
555 8947 : keep_itypat_=.false.
556 8947 : if (present(keep_itypat)) keep_itypat_=keep_itypat
557 8947 : keep_nspden_=.true.
558 8947 : if (present(keep_nspden)) keep_nspden_=keep_nspden
559 :
560 : !Set up parallelism over atoms
561 8947 : paral_atom=(present(comm_atom));if (paral_atom) paral_atom=(xmpi_comm_size(comm_atom)>1)
562 8947 : nullify(my_atmtab);if (present(mpi_atmtab)) my_atmtab => mpi_atmtab
563 8947 : my_comm_atom=xmpi_comm_self;if (present(comm_atom)) my_comm_atom=comm_atom
564 8947 : my_atmtab_allocated=.false.
565 :
566 : !Determine in which case we are (parallelism, ...)
567 : !No parallelism: a single copy operation
568 8947 : paral_case=0;nrhoij_max=nrhoij_in
569 8947 : pawrhoij_out => pawrhoij_cpy
570 8947 : if (paral_atom) then
571 1770 : if (nrhoij_out<nrhoij_in) then ! Parallelism: the copy operation is a scatter
572 56 : call get_my_natom(my_comm_atom,my_nrhoij,nrhoij_in)
573 56 : if (my_nrhoij==nrhoij_out) then
574 56 : call get_my_atmtab(my_comm_atom,my_atmtab,my_atmtab_allocated,paral_atom,nrhoij_in)
575 56 : paral_case=1;nrhoij_max=nrhoij_out
576 56 : pawrhoij_out => pawrhoij_cpy
577 : else
578 0 : msg=' nrhoij_out should be equal to my_natom !'
579 0 : LIBPAW_BUG(msg)
580 : end if
581 : else ! Parallelism: the copy operation is a gather
582 1714 : call get_my_natom(my_comm_atom,my_nrhoij,nrhoij_out)
583 1714 : if (my_nrhoij==nrhoij_in) then
584 1714 : paral_case=2;nrhoij_max=nrhoij_in
585 6799 : LIBPAW_DATATYPE_ALLOCATE(pawrhoij_out,(nrhoij_in))
586 1714 : call pawrhoij_nullify(pawrhoij_out)
587 1714 : if (nrhoij_in>0) then
588 3864 : LIBPAW_ALLOCATE(typat,(nrhoij_in))
589 2576 : LIBPAW_ALLOCATE(nlmn,(nrhoij_in))
590 2945 : do irhoij=1,nrhoij_in
591 2945 : typat(irhoij)=irhoij;nlmn(irhoij)=pawrhoij_in(irhoij)%lmn_size
592 : end do
593 : call pawrhoij_alloc(pawrhoij_out,pawrhoij_cpy(1)%cplex_rhoij,&
594 : & pawrhoij_cpy(1)%nspden,pawrhoij_cpy(1)%nspinor,pawrhoij_cpy(1)%nsppol,typat,&
595 : & lmnsize=nlmn,ngrhoij=pawrhoij_cpy(1)%ngrhoij,nlmnmix=pawrhoij_cpy(1)%lmnmix_sz,&
596 : & qphase=pawrhoij_cpy(1)%qphase,&
597 : & use_rhoij_=pawrhoij_cpy(1)%use_rhoij_,&
598 : & use_rhoijp=pawrhoij_cpy(1)%use_rhoijp,&
599 1288 : & use_rhoijres=pawrhoij_cpy(1)%use_rhoijres)
600 1288 : LIBPAW_DEALLOCATE(typat)
601 1288 : LIBPAW_DEALLOCATE(nlmn)
602 : end if
603 : else
604 0 : msg=' nrhoij_in should be equal to my_natom!'
605 0 : LIBPAW_BUG(msg)
606 : end if
607 : end if
608 : end if
609 :
610 : !Loop on rhoij components
611 8947 : if (nrhoij_max>0) then
612 28896 : do irhoij=1,nrhoij_max
613 20403 : jrhoij=irhoij;if (paral_case==1) jrhoij=my_atmtab(irhoij)
614 :
615 20403 : lmn2_size_in=pawrhoij_in(jrhoij)%lmn2_size
616 20403 : lmn2_size_out=lmn2_size_in
617 20403 : cplex_in=pawrhoij_in(jrhoij)%cplex_rhoij
618 20403 : cplex_out=cplex_in;if(keep_cplex_)cplex_out=pawrhoij_out(irhoij)%cplex_rhoij
619 20403 : qphase_in=pawrhoij_in(jrhoij)%qphase
620 20403 : qphase_out=qphase_in;if(keep_qphase_)qphase_out=pawrhoij_out(irhoij)%qphase
621 20403 : nspden_in=pawrhoij_in(jrhoij)%nspden
622 20403 : nselect=pawrhoij_in(jrhoij)%nrhoijsel
623 20403 : nselect_out=pawrhoij_out(irhoij)%nrhoijsel
624 20403 : nspden_out=nspden_in;if(keep_nspden_)nspden_out=pawrhoij_out(irhoij)%nspden
625 :
626 : change_dim=(pawrhoij_out(irhoij)%cplex_rhoij/=cplex_out.or. &
627 : & pawrhoij_out(irhoij)%qphase/=qphase_out.or. &
628 : & pawrhoij_out(irhoij)%lmn2_size/=lmn2_size_out.or. &
629 : & pawrhoij_out(irhoij)%nspden/=nspden_out.or. &
630 20403 : & nselect/=nselect_out)
631 20403 : cplex=min(cplex_in,cplex_out)
632 20403 : qphase=min(qphase_in,qphase_out)
633 :
634 : ! Scalars
635 20403 : pawrhoij_out(irhoij)%cplex_rhoij=cplex_out+0
636 20403 : pawrhoij_out(irhoij)%qphase=qphase_out+0
637 20403 : pawrhoij_out(irhoij)%nspden=nspden_out+0
638 20403 : pawrhoij_out(irhoij)%lmn2_size=lmn2_size_out+0
639 20403 : pawrhoij_out(irhoij)%lmn_size=pawrhoij_in(jrhoij)%lmn_size+0
640 20403 : if(.not.keep_itypat_) pawrhoij_out(irhoij)%itypat =pawrhoij_in(jrhoij)%itypat+0
641 20403 : if(.not.keep_nspden_) pawrhoij_out(irhoij)%nsppol =pawrhoij_in(jrhoij)%nsppol+0
642 48 : if(.not.keep_nspden_) pawrhoij_out(irhoij)%nspinor=pawrhoij_in(jrhoij)%nspinor+0
643 20403 : pawrhoij_out(irhoij)%nrhoijsel=nselect+0
644 : ! if (pawrhoij_out(irhoij)%itypat/=pawrhoij_in(jrhoij)%itypat) then
645 : ! write(unit=msg,fmt='(a,i3,a)') 'Type of atom ',jrhoij,' is different (dont copy it) !'
646 : ! LIBPAW_COMMENT(msg)
647 : ! end if
648 :
649 : ! Optional pointer: non-zero elements of rhoij
650 20403 : use_rhoijp=pawrhoij_in(jrhoij)%use_rhoijp
651 20403 : if (pawrhoij_out(irhoij)%use_rhoijp/=use_rhoijp) then
652 48 : if (pawrhoij_out(irhoij)%use_rhoijp>0) then
653 0 : LIBPAW_DEALLOCATE(pawrhoij_out(irhoij)%rhoijp)
654 0 : LIBPAW_DEALLOCATE(pawrhoij_out(irhoij)%rhoijselect)
655 : end if
656 48 : if (use_rhoijp>0) then
657 192 : LIBPAW_ALLOCATE(pawrhoij_out(irhoij)%rhoijp,(cplex_out*qphase_out*lmn2_size_out,nspden_out))
658 144 : LIBPAW_ALLOCATE(pawrhoij_out(irhoij)%rhoijselect,(lmn2_size_out))
659 1776 : pawrhoij_out(irhoij)%rhoijselect=0
660 : end if
661 : end if
662 20403 : pawrhoij_out(irhoij)%use_rhoijp=use_rhoijp
663 20403 : if (use_rhoijp>0) then
664 20403 : if (change_dim) then
665 10115 : if(allocated(pawrhoij_out(irhoij)%rhoijp)) then
666 10115 : LIBPAW_DEALLOCATE(pawrhoij_out(irhoij)%rhoijp)
667 : end if
668 40460 : LIBPAW_ALLOCATE(pawrhoij_out(irhoij)%rhoijp,(cplex_out*qphase_out*lmn2_size_out,nspden_out))
669 : end if
670 1988703 : pawrhoij_out(irhoij)%rhoijp(:,:)=zero
671 20403 : if (nspden_out==1) then
672 15470 : if (nspden_in==2) then
673 0 : do iphase=1,qphase
674 0 : i_in=(iphase-1)*lmn2_size_in;i_out=(iphase-1)*lmn2_size_out
675 0 : do ilmn=1,nselect
676 : pawrhoij_out(irhoij)%rhoijp(i_out+1:i_out+cplex,1)= &
677 : & pawrhoij_in(jrhoij)%rhoijp(i_in+1:i_in+cplex,1) &
678 0 : & +pawrhoij_in(jrhoij)%rhoijp(i_in+1:i_in+cplex,2)+zero
679 0 : i_in=i_in+cplex_in;i_out=i_out+cplex_out
680 : end do
681 : end do
682 : else ! nspden_in==1 or 4
683 31020 : do iphase=1,qphase
684 15550 : i_in=(iphase-1)*lmn2_size_in;i_out=(iphase-1)*lmn2_size_out
685 346781 : do ilmn=1,nselect
686 : pawrhoij_out(irhoij)%rhoijp(i_out+1:i_out+cplex,1)= &
687 634478 : & pawrhoij_in(jrhoij)%rhoijp(i_in+1:i_in+cplex,1)+zero
688 331311 : i_in=i_in+cplex_in;i_out=i_out+cplex_out
689 : end do
690 : end do
691 : end if
692 4933 : else if (nspden_out==2) then
693 4158 : if (nspden_in==1) then
694 0 : do iphase=1,qphase
695 0 : i_in=(iphase-1)*lmn2_size_in;i_out=(iphase-1)*lmn2_size_out
696 0 : do ilmn=1,nselect
697 : pawrhoij_out(irhoij)%rhoijp(i_out+1:i_out+cplex,1)= &
698 0 : & half*pawrhoij_in(jrhoij)%rhoijp(i_in+1:i_in+cplex,1)+zero
699 : pawrhoij_out(irhoij)%rhoijp(i_out+1:i_out+cplex,2)= &
700 0 : & half*pawrhoij_in(jrhoij)%rhoijp(i_in+1:i_in+cplex,1)+zero
701 0 : i_in=i_in+cplex_in;i_out=i_out+cplex_out
702 : end do
703 : end do
704 4158 : else if (nspden_in==2) then
705 12474 : do ispden=1,nspden_out
706 20806 : do iphase=1,qphase
707 8332 : i_in=(iphase-1)*lmn2_size_in;i_out=(iphase-1)*lmn2_size_out
708 302100 : do ilmn=1,nselect
709 : pawrhoij_out(irhoij)%rhoijp(i_out+1:i_out+cplex,ispden)= &
710 571664 : & pawrhoij_in(jrhoij)%rhoijp(i_in+1:i_in+cplex,ispden)+zero
711 293784 : i_in=i_in+cplex_in;i_out=i_out+cplex_out
712 : end do
713 : end do
714 : end do
715 : else ! nspden_in==4
716 0 : do iphase=1,qphase
717 0 : i_in=(iphase-1)*lmn2_size_in;i_out=(iphase-1)*lmn2_size_out
718 0 : do ilmn=1,nselect
719 : pawrhoij_out(irhoij)%rhoijp(i_out+1:i_out+cplex,1)= &
720 : & half*(pawrhoij_in(jrhoij)%rhoijp(i_in+1:i_in+cplex,1) &
721 0 : & +pawrhoij_in(jrhoij)%rhoijp(i_in+1:i_in+cplex,4))+zero
722 : pawrhoij_out(irhoij)%rhoijp(i_out+1:i_out+cplex,2)= &
723 : & half*(pawrhoij_in(jrhoij)%rhoijp(i_in+1:i_in+cplex,1) &
724 0 : & -pawrhoij_in(jrhoij)%rhoijp(i_in+1:i_in+cplex,4))+zero
725 0 : i_in=i_in+cplex_in;i_out=i_out+cplex_out
726 : end do
727 : end do
728 : end if
729 775 : else if (nspden_out==4) then
730 775 : if (nspden_in==1) then
731 16 : do iphase=1,qphase
732 8 : i_in=(iphase-1)*lmn2_size_in;i_out=(iphase-1)*lmn2_size_out
733 520 : do ilmn=1,nselect
734 : pawrhoij_out(irhoij)%rhoijp(i_out+1:i_out+cplex,1)= &
735 1008 : & pawrhoij_in(jrhoij)%rhoijp(i_in+1:i_in+cplex,1)+zero
736 3528 : pawrhoij_out(irhoij)%rhoijp(i_out+1:i_out+cplex,2:4)=zero
737 512 : i_in=i_in+cplex_in;i_out=i_out+cplex_out
738 : end do
739 : end do
740 767 : else if (nspden_in==2) then
741 2 : do iphase=1,qphase
742 1 : i_in=(iphase-1)*lmn2_size_in;i_out=(iphase-1)*lmn2_size_out
743 29 : do ilmn=1,nselect
744 : pawrhoij_out(irhoij)%rhoijp(i_out+1:i_out+cplex,1)= &
745 : & pawrhoij_in(jrhoij)%rhoijp(i_in+1:i_in+cplex,1) &
746 54 : & +pawrhoij_in(jrhoij)%rhoijp(i_in+1:i_in+cplex,2)+zero
747 : pawrhoij_out(irhoij)%rhoijp(i_out+1:i_out+cplex,4)= &
748 : & pawrhoij_in(jrhoij)%rhoijp(i_in+1:i_in+cplex,1) &
749 54 : & -pawrhoij_in(jrhoij)%rhoijp(i_in+1:i_in+cplex,2)+zero
750 135 : pawrhoij_out(irhoij)%rhoijp(i_out+1:i_out+cplex,2:3)=zero
751 28 : i_in=i_in+cplex_in;i_out=i_out+cplex_out
752 : end do
753 : end do
754 : else ! nspden_in==4
755 3830 : do ispden=1,nspden_out
756 6894 : do iphase=1,qphase
757 3064 : i_in=(iphase-1)*lmn2_size_in;i_out=(iphase-1)*lmn2_size_out
758 191352 : do ilmn=1,nselect
759 : pawrhoij_out(irhoij)%rhoijp(i_out+1:i_out+cplex,ispden)= &
760 527468 : & pawrhoij_in(jrhoij)%rhoijp(i_in+1:i_in+cplex,ispden)+zero
761 188288 : i_in=i_in+cplex_in;i_out=i_out+cplex_out
762 : end do
763 : end do
764 : end do
765 : end if
766 : end if
767 : end if
768 :
769 : ! Optional pointer: indexes for non-zero elements selection
770 : if (use_rhoijp>0) then
771 20403 : if (change_dim) then
772 10115 : if(allocated(pawrhoij_out(irhoij)%rhoijselect)) then
773 10115 : LIBPAW_DEALLOCATE(pawrhoij_out(irhoij)%rhoijselect)
774 : end if
775 30345 : LIBPAW_ALLOCATE(pawrhoij_out(irhoij)%rhoijselect,(lmn2_size_out))
776 : end if
777 1106649 : pawrhoij_out(irhoij)%rhoijselect=0
778 521959 : pawrhoij_out(irhoij)%rhoijselect(1:nselect)=pawrhoij_in(jrhoij)%rhoijselect(1:nselect)+0
779 : end if
780 :
781 : ! Optional pointer: indexes of rhoij to be mixed
782 20403 : lmnmix=pawrhoij_in(jrhoij)%lmnmix_sz
783 20403 : if (pawrhoij_out(irhoij)%lmnmix_sz/=lmnmix) then
784 5036 : if (pawrhoij_out(irhoij)%lmnmix_sz>0) then
785 2493 : LIBPAW_DEALLOCATE(pawrhoij_out(irhoij)%kpawmix)
786 : end if
787 5036 : if (lmnmix>0) then
788 7629 : LIBPAW_ALLOCATE(pawrhoij_out(irhoij)%kpawmix,(lmnmix))
789 164911 : pawrhoij_out(irhoij)%kpawmix=0
790 : end if
791 5036 : pawrhoij_out(irhoij)%lmnmix_sz=lmnmix
792 : end if
793 302556 : if (lmnmix>0) pawrhoij_out(irhoij)%kpawmix(1:lmnmix)=pawrhoij_in(jrhoij)%kpawmix(1:lmnmix)
794 :
795 : ! Optional pointer: gradients of rhoij
796 20403 : ngrhoij=pawrhoij_in(jrhoij)%ngrhoij
797 20403 : if (pawrhoij_out(irhoij)%ngrhoij/=ngrhoij) then
798 0 : if (pawrhoij_out(irhoij)%ngrhoij>0) then
799 0 : LIBPAW_DEALLOCATE(pawrhoij_out(irhoij)%grhoij)
800 : end if
801 0 : if (ngrhoij>0) then
802 0 : LIBPAW_ALLOCATE(pawrhoij_out(irhoij)%grhoij,(ngrhoij,cplex_out*qphase_out*lmn2_size_out,nspden_out))
803 : end if
804 0 : pawrhoij_out(irhoij)%ngrhoij=ngrhoij
805 : end if
806 20403 : if (ngrhoij>0) then
807 0 : if (change_dim) then
808 0 : LIBPAW_DEALLOCATE(pawrhoij_out(irhoij)%grhoij)
809 0 : LIBPAW_ALLOCATE(pawrhoij_out(irhoij)%grhoij,(ngrhoij,cplex_out*qphase_out*lmn2_size_out,nspden_out))
810 : end if
811 0 : pawrhoij_out(irhoij)%grhoij(:,:,:)=zero
812 0 : if (nspden_out==1) then
813 0 : if (nspden_in==2) then
814 0 : do iphase=1,qphase
815 0 : i_in=(iphase-1)*lmn2_size_in;i_out=(iphase-1)*lmn2_size_out
816 0 : do ilmn=1,lmn2_size_out
817 : pawrhoij_out(irhoij)%grhoij(1:ngrhoij,i_out+1:i_out+cplex,1)= &
818 : & pawrhoij_in(jrhoij)%grhoij(1:ngrhoij,i_in+1:i_in+cplex,1) &
819 0 : & +pawrhoij_in(jrhoij)%grhoij(1:ngrhoij,i_in+1:i_in+cplex,2)+zero
820 0 : i_in=i_in+cplex_in;i_out=i_out+cplex_out
821 : end do
822 : end do
823 : else ! nspden_in==1 or 4
824 0 : do iphase=1,qphase
825 0 : i_in=(iphase-1)*lmn2_size_in;i_out=(iphase-1)*lmn2_size_out
826 0 : do ilmn=1,lmn2_size_out
827 : pawrhoij_out(irhoij)%grhoij(1:ngrhoij,i_out+1:i_out+cplex,1)= &
828 0 : & pawrhoij_in(jrhoij)%grhoij(1:ngrhoij,i_in+1:i_in+cplex,1)+zero
829 0 : i_in=i_in+cplex_in;i_out=i_out+cplex_out
830 : end do
831 : end do
832 : end if
833 0 : else if (nspden_out==2) then
834 0 : if (nspden_in==1) then
835 0 : do iphase=1,qphase
836 0 : i_in=(iphase-1)*lmn2_size_in;i_out=(iphase-1)*lmn2_size_out
837 0 : do ilmn=1,lmn2_size_out
838 : pawrhoij_out(irhoij)%grhoij(1:ngrhoij,i_out+1:i_out+cplex,1)= &
839 0 : & half*pawrhoij_in(jrhoij)%grhoij(1:ngrhoij,i_in+1:i_in+cplex,1)+zero
840 : pawrhoij_out(irhoij)%grhoij(1:ngrhoij,i_out+1:i_out+cplex,2)= &
841 0 : & half*pawrhoij_in(jrhoij)%grhoij(1:ngrhoij,i_in+1:i_in+cplex,1)+zero
842 0 : i_in=i_in+cplex_in;i_out=i_out+cplex_out
843 : end do
844 : end do
845 0 : else if (nspden_in==2) then
846 0 : do ispden=1,nspden_out
847 0 : do iphase=1,qphase
848 0 : i_in=(iphase-1)*lmn2_size_in;i_out=(iphase-1)*lmn2_size_out
849 0 : do ilmn=1,lmn2_size_out
850 : pawrhoij_out(irhoij)%grhoij(1:ngrhoij,i_out+1:i_out+cplex,ispden)= &
851 0 : & pawrhoij_in(jrhoij)%grhoij(1:ngrhoij,i_in+1:i_in+cplex,ispden)+zero
852 0 : i_in=i_in+cplex_in;i_out=i_out+cplex_out
853 : end do
854 : end do
855 : end do
856 : else ! nspden_in==4
857 0 : do iphase=1,qphase
858 0 : i_in=(iphase-1)*lmn2_size_in;i_out=(iphase-1)*lmn2_size_out
859 0 : do ilmn=1,lmn2_size_out
860 : pawrhoij_out(irhoij)%grhoij(1:ngrhoij,i_out+1:i_out+cplex,1)= &
861 : & half*(pawrhoij_in(jrhoij)%grhoij(1:ngrhoij,i_in+1:i_in+cplex,1) &
862 0 : & +pawrhoij_in(jrhoij)%grhoij(1:ngrhoij,i_in+1:i_in+cplex,4))+zero
863 : pawrhoij_out(irhoij)%grhoij(1:ngrhoij,i_out+1:i_out+cplex,2)= &
864 : & half*(pawrhoij_in(jrhoij)%grhoij(1:ngrhoij,i_in+1:i_in+cplex,1) &
865 0 : & -pawrhoij_in(jrhoij)%grhoij(1:ngrhoij,i_in+1:i_in+cplex,4))+zero
866 0 : i_in=i_in+cplex_in;i_out=i_out+cplex_out
867 : end do
868 : end do
869 : end if
870 0 : else if (nspden_out==4) then
871 0 : if (nspden_in==1) then
872 0 : do iphase=1,qphase
873 0 : i_in=(iphase-1)*lmn2_size_in;i_out=(iphase-1)*lmn2_size_out
874 0 : do ilmn=1,lmn2_size_out
875 : pawrhoij_out(irhoij)%grhoij(1:ngrhoij,i_out+1:i_out+cplex,1)= &
876 0 : & pawrhoij_in(jrhoij)%grhoij(1:ngrhoij,i_in+1:i_in+cplex,1)+zero
877 0 : pawrhoij_out(irhoij)%grhoij(1:ngrhoij,i_out+1:i_out+cplex,2:4)=zero
878 0 : i_in=i_in+cplex_in;i_out=i_out+cplex_out
879 : end do
880 : end do
881 0 : else if (nspden_in==2) then
882 0 : do iphase=1,qphase
883 0 : i_in=(iphase-1)*lmn2_size_in;i_out=(iphase-1)*lmn2_size_out
884 0 : do ilmn=1,lmn2_size_out
885 : pawrhoij_out(irhoij)%grhoij(1:ngrhoij,i_out+1:i_out+cplex,1)= &
886 : & pawrhoij_in(jrhoij)%grhoij(1:ngrhoij,i_in+1:i_in+cplex,1) &
887 0 : & +pawrhoij_in(jrhoij)%grhoij(1:ngrhoij,i_in+1:i_in+cplex,2)+zero
888 : pawrhoij_out(irhoij)%grhoij(1:ngrhoij,i_out+1:i_out+cplex,4)= &
889 : & pawrhoij_in(jrhoij)%grhoij(1:ngrhoij,i_in+1:i_in+cplex,1) &
890 0 : & -pawrhoij_in(jrhoij)%grhoij(1:ngrhoij,i_in+1:i_in+cplex,2)+zero
891 0 : pawrhoij_out(irhoij)%grhoij(1:ngrhoij,i_out+1:i_out+cplex,2:3)=zero
892 0 : i_in=i_in+cplex_in;i_out=i_out+cplex_out
893 : end do
894 : end do
895 : else ! nspden_in==4
896 0 : do ispden=1,nspden_out
897 0 : do iphase=1,qphase
898 0 : i_in=(iphase-1)*lmn2_size_in;i_out=(iphase-1)*lmn2_size_out
899 0 : do ilmn=1,lmn2_size_out
900 : pawrhoij_out(irhoij)%grhoij(1:ngrhoij,i_out+1:i_out+cplex,ispden)= &
901 0 : & pawrhoij_in(jrhoij)%grhoij(1:ngrhoij,i_in+1:i_in+cplex,ispden)+zero
902 0 : i_in=i_in+cplex_in;i_out=i_out+cplex_out
903 : end do
904 : end do
905 : end do
906 : end if
907 : end if
908 : end if
909 :
910 : ! Optional pointer: residuals of rhoij
911 20403 : use_rhoijres=pawrhoij_in(jrhoij)%use_rhoijres
912 20403 : if (pawrhoij_out(irhoij)%use_rhoijres/=use_rhoijres) then
913 5036 : if (pawrhoij_out(irhoij)%use_rhoijres>0) then
914 2493 : LIBPAW_DEALLOCATE(pawrhoij_out(irhoij)%rhoijres)
915 : end if
916 5036 : if (use_rhoijres>0) then
917 10172 : LIBPAW_ALLOCATE(pawrhoij_out(irhoij)%rhoijres,(cplex_out*qphase_out*lmn2_size_out,nspden_out))
918 : end if
919 5036 : pawrhoij_out(irhoij)%use_rhoijres=use_rhoijres
920 : end if
921 20403 : if (use_rhoijres>0) then
922 5686 : if (change_dim) then
923 2547 : LIBPAW_DEALLOCATE(pawrhoij_out(irhoij)%rhoijres)
924 10188 : LIBPAW_ALLOCATE(pawrhoij_out(irhoij)%rhoijres,(cplex_out*qphase_out*lmn2_size_out,nspden_out))
925 : end if
926 480189 : pawrhoij_out(irhoij)%rhoijres(:,:)=zero
927 5686 : if (nspden_out==1) then
928 5030 : if (nspden_in==2) then
929 0 : do iphase=1,qphase
930 0 : i_in=(iphase-1)*lmn2_size_in;i_out=(iphase-1)*lmn2_size_out
931 0 : do ilmn=1,lmn2_size_out
932 : pawrhoij_out(irhoij)%rhoijres(i_out+1:i_out+cplex,1)= &
933 : & pawrhoij_in(jrhoij)%rhoijres(i_in+1:i_in+cplex,1) &
934 0 : & +pawrhoij_in(jrhoij)%rhoijres(i_in+1:i_in+cplex,2)+zero
935 0 : i_in=i_in+cplex_in;i_out=i_out+cplex_out
936 : end do
937 : end do
938 : else ! nspden_in==1 or 4
939 10060 : do iphase=1,qphase
940 5030 : i_in=(iphase-1)*lmn2_size_in;i_out=(iphase-1)*lmn2_size_out
941 214307 : do ilmn=1,lmn2_size_out
942 : pawrhoij_out(irhoij)%rhoijres(i_out+1:i_out+cplex,1)= &
943 409286 : & pawrhoij_in(jrhoij)%rhoijres(i_in+1:i_in+cplex,1)+zero
944 209277 : i_in=i_in+cplex_in;i_out=i_out+cplex_out
945 : end do
946 : end do
947 : end if
948 656 : else if (nspden_out==2) then
949 476 : if (nspden_in==1) then
950 0 : do iphase=1,qphase
951 0 : i_in=(iphase-1)*lmn2_size_in;i_out=(iphase-1)*lmn2_size_out
952 0 : do ilmn=1,lmn2_size_out
953 : pawrhoij_out(irhoij)%rhoijres(i_out+1:i_out+cplex,1)= &
954 0 : & half*pawrhoij_in(jrhoij)%rhoijres(i_in+1:i_in+cplex,1)+zero
955 : pawrhoij_out(irhoij)%rhoijres(i_out+1:i_out+cplex,2)= &
956 0 : & half*pawrhoij_in(jrhoij)%rhoijres(i_in+1:i_in+cplex,1)+zero
957 0 : i_in=i_in+cplex_in;i_out=i_out+cplex_out
958 : end do
959 : end do
960 476 : else if (nspden_in==2) then
961 1428 : do ispden=1,nspden_out
962 2380 : do iphase=1,qphase
963 952 : i_in=(iphase-1)*lmn2_size_in;i_out=(iphase-1)*lmn2_size_out
964 117230 : do ilmn=1,lmn2_size_out
965 : pawrhoij_out(irhoij)%rhoijres(i_out+1:i_out+cplex,ispden)= &
966 230796 : & pawrhoij_in(jrhoij)%rhoijres(i_in+1:i_in+cplex,ispden)+zero
967 116278 : i_in=i_in+cplex_in;i_out=i_out+cplex_out
968 : end do
969 : end do
970 : end do
971 : else ! nspden_in==4
972 0 : do iphase=1,qphase
973 0 : i_in=(iphase-1)*lmn2_size_in;i_out=(iphase-1)*lmn2_size_out
974 0 : do ilmn=1,lmn2_size_out
975 : pawrhoij_out(irhoij)%rhoijres(i_out+1:i_out+cplex,1)= &
976 : & half*(pawrhoij_in(jrhoij)%rhoijres(i_in+1:i_in+cplex,1) &
977 0 : & +pawrhoij_in(jrhoij)%rhoijres(i_in+1:i_in+cplex,4))+zero
978 : pawrhoij_out(irhoij)%rhoijres(i_out+1:i_out+cplex,2)= &
979 : & half*(pawrhoij_in(jrhoij)%rhoijres(i_in+1:i_in+cplex,1) &
980 0 : & -pawrhoij_in(jrhoij)%rhoijres(i_in+1:i_in+cplex,4))+zero
981 0 : i_in=i_in+cplex_in;i_out=i_out+cplex_out
982 : end do
983 : end do
984 : end if
985 180 : else if (nspden_out==4) then
986 180 : if (nspden_in==1) then
987 0 : do iphase=1,qphase
988 0 : i_in=(iphase-1)*lmn2_size_in;i_out=(iphase-1)*lmn2_size_out
989 0 : do ilmn=1,lmn2_size_out
990 : pawrhoij_out(irhoij)%rhoijres(i_out+1:i_out+cplex,1)= &
991 0 : & pawrhoij_in(jrhoij)%rhoijres(i_in+1:i_in+cplex,1)+zero
992 0 : pawrhoij_out(irhoij)%rhoijres(i_out+1:i_out+cplex,2:4)=zero
993 0 : i_in=i_in+cplex_in;i_out=i_out+cplex_out
994 : end do
995 : end do
996 180 : else if (nspden_in==2) then
997 0 : do iphase=1,qphase
998 0 : i_in=(iphase-1)*lmn2_size_in;i_out=(iphase-1)*lmn2_size_out
999 0 : do ilmn=1,lmn2_size_out
1000 : pawrhoij_out(irhoij)%rhoijres(i_out+1:i_out+cplex,1)= &
1001 : & pawrhoij_in(jrhoij)%rhoijres(i_in+1:i_in+cplex,1) &
1002 0 : & +pawrhoij_in(jrhoij)%rhoijres(i_in+1:i_in+cplex,2)+zero
1003 : pawrhoij_out(irhoij)%rhoijres(i_out+1:i_out+cplex,4)= &
1004 : & pawrhoij_in(jrhoij)%rhoijres(i_in+1:i_in+cplex,1) &
1005 0 : & -pawrhoij_in(jrhoij)%rhoijres(i_in+1:i_in+cplex,2)+zero
1006 0 : pawrhoij_out(irhoij)%rhoijres(i_out+1:i_out+cplex,2:3)=zero
1007 0 : i_in=i_in+cplex_in;i_out=i_out+cplex_out
1008 : end do
1009 : end do
1010 : else ! nspden_in==4
1011 900 : do ispden=1,nspden_out
1012 1620 : do iphase=1,qphase
1013 720 : i_in=(iphase-1)*lmn2_size_in;i_out=(iphase-1)*lmn2_size_out
1014 82952 : do ilmn=1,lmn2_size_out
1015 : pawrhoij_out(irhoij)%rhoijres(i_out+1:i_out+cplex,ispden)= &
1016 228804 : & pawrhoij_in(jrhoij)%rhoijres(i_in+1:i_in+cplex,ispden)+zero
1017 82232 : i_in=i_in+cplex_in;i_out=i_out+cplex_out
1018 : end do
1019 : end do
1020 : end do
1021 : end if
1022 : end if
1023 : end if
1024 :
1025 : ! Optional pointer: non-symmetrized rhoij
1026 20403 : use_rhoij_=pawrhoij_in(jrhoij)%use_rhoij_
1027 20403 : if (pawrhoij_out(irhoij)%use_rhoij_/=use_rhoij_) then
1028 0 : if (pawrhoij_out(irhoij)%use_rhoij_>0) then
1029 0 : LIBPAW_DEALLOCATE(pawrhoij_out(irhoij)%rhoij_)
1030 : end if
1031 0 : if (use_rhoij_>0) then
1032 0 : LIBPAW_ALLOCATE(pawrhoij_out(irhoij)%rhoij_,(cplex_out*qphase_out*lmn2_size_out,nspden_out))
1033 : end if
1034 0 : pawrhoij_out(irhoij)%use_rhoij_=use_rhoij_
1035 : end if
1036 28896 : if (use_rhoij_>0) then
1037 0 : if (change_dim) then
1038 0 : if(allocated(pawrhoij_out(irhoij)%rhoij_)) then
1039 0 : LIBPAW_DEALLOCATE(pawrhoij_out(irhoij)%rhoij_)
1040 : end if
1041 0 : LIBPAW_ALLOCATE(pawrhoij_out(irhoij)%rhoij_,(cplex_out*qphase_out*lmn2_size_out,nspden_out))
1042 : end if
1043 0 : pawrhoij_out(irhoij)%rhoij_(:,:)=zero
1044 0 : if (nspden_out==1) then
1045 0 : if (nspden_in==2) then
1046 0 : do iphase=1,qphase
1047 0 : i_in=(iphase-1)*lmn2_size_in;i_out=(iphase-1)*lmn2_size_out
1048 0 : do ilmn=1,lmn2_size_out
1049 : pawrhoij_out(irhoij)%rhoij_(i_out+1:i_out+cplex,1)= &
1050 : & pawrhoij_in(jrhoij)%rhoij_(i_in+1:i_in+cplex,1) &
1051 0 : & +pawrhoij_in(jrhoij)%rhoij_(i_in+1:i_in+cplex,2)+zero
1052 0 : i_in=i_in+cplex_in;i_out=i_out+cplex_out
1053 : end do
1054 : end do
1055 : else ! nspden_in==1 or 4
1056 0 : do iphase=1,qphase
1057 0 : i_in=(iphase-1)*lmn2_size_in;i_out=(iphase-1)*lmn2_size_out
1058 0 : do ilmn=1,lmn2_size_out
1059 : pawrhoij_out(irhoij)%rhoij_(i_out+1:i_out+cplex,1)= &
1060 0 : & pawrhoij_in(jrhoij)%rhoij_(i_in+1:i_in+cplex,1)+zero
1061 0 : i_in=i_in+cplex_in;i_out=i_out+cplex_out
1062 : end do
1063 : end do
1064 : end if
1065 0 : else if (nspden_out==2) then
1066 0 : if (nspden_in==1) then
1067 0 : do iphase=1,qphase
1068 0 : i_in=(iphase-1)*lmn2_size_in;i_out=(iphase-1)*lmn2_size_out
1069 0 : do ilmn=1,lmn2_size_out
1070 : pawrhoij_out(irhoij)%rhoij_(i_out+1:i_out+cplex,1)= &
1071 0 : & half*pawrhoij_in(jrhoij)%rhoij_(i_in+1:i_in+cplex,1)+zero
1072 : pawrhoij_out(irhoij)%rhoij_(i_out+1:i_out+cplex,2)= &
1073 0 : & half*pawrhoij_in(irhoij)%rhoij_(i_in+1:i_in+cplex,1)+zero
1074 0 : i_in=i_in+cplex_in;i_out=i_out+cplex_out
1075 : end do
1076 : end do
1077 0 : else if (nspden_in==2) then
1078 0 : do ispden=1,nspden_out
1079 0 : do iphase=1,qphase
1080 0 : i_in=(iphase-1)*lmn2_size_in;i_out=(iphase-1)*lmn2_size_out
1081 0 : do ilmn=1,lmn2_size_out
1082 : pawrhoij_out(irhoij)%rhoij_(i_out+1:i_out+cplex,ispden)= &
1083 0 : & pawrhoij_in(jrhoij)%rhoij_(i_in+1:i_in+cplex,ispden)+zero
1084 0 : i_in=i_in+cplex_in;i_out=i_out+cplex_out
1085 : end do
1086 : end do
1087 : end do
1088 : else ! nspden_in==4
1089 0 : do iphase=1,qphase
1090 0 : i_in=(iphase-1)*lmn2_size_in;i_out=(iphase-1)*lmn2_size_out
1091 0 : do ilmn=1,lmn2_size_out
1092 : pawrhoij_out(irhoij)%rhoij_(i_out+1:i_out+cplex,1)= &
1093 : & half*(pawrhoij_in(jrhoij)%rhoij_(i_in+1:i_in+cplex,1) &
1094 0 : & +pawrhoij_in(jrhoij)%rhoij_(i_in+1:i_in+cplex,4))+zero
1095 : pawrhoij_out(irhoij)%rhoij_(i_out+1:i_out+cplex,2)= &
1096 : & half*(pawrhoij_in(jrhoij)%rhoij_(i_in+1:i_in+cplex,1) &
1097 0 : & -pawrhoij_in(jrhoij)%rhoij_(i_in+1:i_in+cplex,4))+zero
1098 0 : i_in=i_in+cplex_in;i_out=i_out+cplex_out
1099 : end do
1100 : end do
1101 : end if
1102 0 : else if (nspden_out==4) then
1103 0 : if (nspden_in==1) then
1104 0 : do iphase=1,qphase
1105 0 : i_in=(iphase-1)*lmn2_size_in;i_out=(iphase-1)*lmn2_size_out
1106 0 : do ilmn=1,lmn2_size_out
1107 : pawrhoij_out(irhoij)%rhoij_(i_out+1:i_out+cplex,1)= &
1108 0 : & pawrhoij_in(jrhoij)%rhoij_(i_in+1:i_in+cplex,1)+zero
1109 0 : pawrhoij_out(irhoij)%rhoij_(i_out+1:i_out+cplex,2:4)=zero
1110 0 : i_in=i_in+cplex_in;i_out=i_out+cplex_out
1111 : end do
1112 : end do
1113 0 : else if (nspden_in==2) then
1114 0 : do iphase=1,qphase
1115 0 : i_in=(iphase-1)*lmn2_size_in;i_out=(iphase-1)*lmn2_size_out
1116 0 : do ilmn=1,lmn2_size_out
1117 : pawrhoij_out(irhoij)%rhoij_(i_out+1:i_out+cplex,1)= &
1118 : & pawrhoij_in(jrhoij)%rhoij_(i_in+1:i_in+cplex,1) &
1119 0 : & +pawrhoij_in(jrhoij)%rhoij_(i_in+1:i_in+cplex,2)+zero
1120 : pawrhoij_out(irhoij)%rhoij_(i_out+1:i_out+cplex,4)= &
1121 : & pawrhoij_in(jrhoij)%rhoij_(i_in+1:i_in+cplex,1) &
1122 0 : & -pawrhoij_in(jrhoij)%rhoij_(i_in+1:i_in+cplex,2)+zero
1123 0 : pawrhoij_out(irhoij)%rhoij_(i_out+1:i_out+cplex,2:3)=zero
1124 0 : i_in=i_in+cplex_in;i_out=i_out+cplex_out
1125 : end do
1126 : end do
1127 : else ! nspden_in==4
1128 0 : do ispden=1,nspden_out
1129 0 : do iphase=1,qphase
1130 0 : i_in=(iphase-1)*lmn2_size_in;i_out=(iphase-1)*lmn2_size_out
1131 0 : do ilmn=1,lmn2_size_out
1132 : pawrhoij_out(irhoij)%rhoij_(i_out+1:i_out+cplex,ispden)= &
1133 0 : & pawrhoij_in(jrhoij)%rhoij_(i_in+1:i_in+cplex,ispden)+zero
1134 0 : i_in=i_in+cplex_in;i_out=i_out+cplex_out
1135 : end do
1136 : end do
1137 : end do
1138 : end if
1139 : end if
1140 : end if
1141 :
1142 : end do ! irhoij
1143 : end if
1144 :
1145 : !Parallel case: do a gather if needed
1146 8947 : if (paral_case==2) then
1147 1714 : call pawrhoij_free(pawrhoij_cpy)
1148 1714 : call pawrhoij_gather(pawrhoij_out,pawrhoij_cpy,-1,my_comm_atom)
1149 1714 : call pawrhoij_free(pawrhoij_out)
1150 1714 : LIBPAW_DATATYPE_DEALLOCATE(pawrhoij_out)
1151 :
1152 : ! Sequential case: fill missing elements
1153 7233 : else if (paral_case==0) then
1154 7177 : if (nrhoij_in<nrhoij_out) then
1155 44 : do irhoij=nrhoij_in+1,nrhoij_out
1156 32 : pawrhoij_cpy(irhoij)%nrhoijsel=0
1157 32 : if (pawrhoij_cpy(irhoij)%use_rhoijp>0) then
1158 1184 : pawrhoij_cpy(irhoij)%rhoijselect=0
1159 1216 : pawrhoij_cpy(irhoij)%rhoijp=zero
1160 : end if
1161 32 : if (pawrhoij_cpy(irhoij)%lmnmix_sz>0) pawrhoij_cpy(irhoij)%kpawmix=0
1162 32 : if (pawrhoij_cpy(irhoij)%ngrhoij>0) pawrhoij_cpy(irhoij)%grhoij=zero
1163 32 : if (pawrhoij_cpy(irhoij)%use_rhoij_>0) pawrhoij_cpy(irhoij)%rhoij_=zero
1164 44 : if (pawrhoij_cpy(irhoij)%use_rhoijres>0) pawrhoij_cpy(irhoij)%rhoijres=zero
1165 : end do
1166 : end if
1167 : end if
1168 :
1169 : !Destroy atom table used for parallelism
1170 8947 : call free_my_atmtab(my_atmtab,my_atmtab_allocated)
1171 :
1172 17894 : end subroutine pawrhoij_copy
1173 : !!***
1174 :
1175 : !----------------------------------------------------------------------
1176 :
1177 : !!****f* m_pawrhoij/pawrhoij_gather
1178 : !! NAME
1179 : !! pawrhoij_gather
1180 : !!
1181 : !! FUNCTION
1182 : !! (All)Gather pawrhoij datastructures
1183 : !!
1184 : !! INPUTS
1185 : !! master=master communicator receiving data ; if -1 do a ALLGATHER
1186 : !! comm_atom= communicator
1187 : !! pawrhoij_in(:)<type(pawrhoij_type)>= input rhoij datastructures on every process
1188 : !! with_grhoij : optional argument (logical, default=.TRUE.)
1189 : !! TRUE if pawrhoij%grhoij field is included in the gather operation
1190 : !! with_lmnmix : optional argument (logical, default=.TRUE.)
1191 : !! TRUE if pawrhoij%lmnmix field is included in the gather operation
1192 : !! with_rhoijp : optional argument (logical, default=.TRUE.)
1193 : !! TRUE if pawrhoij%rhoijp and pawrhoij%rhoijselect fields
1194 : !! are included in the gather operation
1195 : !! with_rhoijres : optional argument (logical, default=.TRUE.)
1196 : !! TRUE if pawrhoij%rhoijres field is included in the gather operation
1197 : !! with_rhoij_ : optional argument (logical, default=.TRUE.)
1198 : !! TRUE if pawrhoij%rhoij_ field is included in the gather operation
1199 : !!
1200 : !! OUTPUT
1201 : !! pawrhoij_gathered(:)<type(pawrhoij_type)>= output rhoij datastructure
1202 : !!
1203 : !! NOTES
1204 : !! The gathered structure are ordered like in sequential mode.
1205 : !!
1206 : !! SOURCE
1207 :
1208 :
1209 2182 : subroutine pawrhoij_gather(pawrhoij_in,pawrhoij_gathered,master,comm_atom, &
1210 : & with_grhoij,with_lmnmix,with_rhoijp,with_rhoijres,with_rhoij_) ! optional arguments
1211 :
1212 : !Arguments ------------------------------------
1213 : !scalars
1214 : integer,intent(in) :: master,comm_atom
1215 : logical,intent(in),optional :: with_grhoij,with_lmnmix,with_rhoijp,with_rhoijres,with_rhoij_
1216 : !arrays
1217 : type(pawrhoij_type),intent(in) :: pawrhoij_in(:)
1218 : type(pawrhoij_type),intent(inout) :: pawrhoij_gathered(:)
1219 : !Local variables-------------------------------
1220 : !scalars
1221 : integer :: buf_dp_size,buf_dp_size_all,buf_int_size,buf_int_size_all
1222 : integer :: cplex,ierr,ii,indx_dp,indx_int,irhoij,isp,jj,jrhoij,lmn2_size,lmnmix,me_atom
1223 : integer :: ngrhoij,nproc_atom,nrhoij_in,nrhoij_in_sum,nrhoij_out,nselect,nspden
1224 : integer :: qphase,rhoij_size2,use_rhoijp,use_rhoijres,use_rhoij_
1225 : logical :: my_atmtab_allocated,paral_atom
1226 : logical :: with_grhoij_,with_lmnmix_,with_rhoijp_,with_rhoijres_,with_rhoij__
1227 : character(len=500) :: msg
1228 : !arrays
1229 : integer :: bufsz(2)
1230 2182 : integer,allocatable :: buf_int(:),buf_int_all(:)
1231 2182 : integer,allocatable :: count_dp(:),count_int(:),count_tot(:),displ_dp(:),displ_int(:)
1232 2182 : integer, pointer :: my_atmtab(:)
1233 2182 : real(dp),allocatable :: buf_dp(:),buf_dp_all(:)
1234 :
1235 : ! *************************************************************************
1236 :
1237 2182 : nrhoij_in=size(pawrhoij_in);nrhoij_out=size(pawrhoij_gathered)
1238 :
1239 4364 : nproc_atom=xmpi_comm_size(comm_atom)
1240 2182 : me_atom=xmpi_comm_rank(comm_atom)
1241 :
1242 2182 : if (nproc_atom==1) then
1243 0 : if (master==-1.or.me_atom==master) then
1244 0 : call pawrhoij_copy(pawrhoij_in,pawrhoij_gathered,.false.,.false.,.false.)
1245 : end if
1246 0 : return
1247 : end if
1248 :
1249 : !Test on sizes
1250 2182 : nrhoij_in_sum=nrhoij_in
1251 2182 : call xmpi_sum(nrhoij_in_sum,comm_atom,ierr)
1252 2182 : if (master==-1) then
1253 1740 : if (nrhoij_out/=nrhoij_in_sum) then
1254 0 : msg='Wrong sizes sum[nrhoij_ij]/=nrhoij_out !'
1255 0 : LIBPAW_BUG(msg)
1256 : end if
1257 : else
1258 442 : if (me_atom==master.and.nrhoij_out/=nrhoij_in_sum) then
1259 0 : msg='(2) pawrhoij_gathered wrongly allocated !'
1260 0 : LIBPAW_BUG(msg)
1261 : end if
1262 : end if
1263 :
1264 : !Optional arguments
1265 2182 : with_grhoij_ =.true.;if (present(with_grhoij)) with_grhoij_ =with_grhoij
1266 2182 : with_lmnmix_ =.true.;if (present(with_lmnmix)) with_lmnmix_ =with_lmnmix
1267 2182 : with_rhoijp_ =.true.;if (present(with_rhoijp)) with_rhoijp_ =with_rhoijp
1268 2182 : with_rhoijres_=.true.;if (present(with_rhoijres))with_rhoijres_=with_rhoijres
1269 2182 : with_rhoij__ =.true.;if (present(with_rhoij_)) with_rhoij__ =with_rhoij_
1270 :
1271 : !Retrieve table of atoms
1272 2182 : paral_atom=.true.;nullify(my_atmtab)
1273 : call get_my_atmtab(comm_atom,my_atmtab,my_atmtab_allocated,paral_atom,nrhoij_in_sum, &
1274 2182 : & my_natom_ref=nrhoij_in)
1275 :
1276 : !Compute sizes of buffers
1277 2182 : buf_int_size=0;buf_dp_size=0
1278 2182 : nselect=0;lmnmix=0;ngrhoij=0;rhoij_size2=0
1279 2182 : use_rhoijp=0;use_rhoijres=0;use_rhoij_=0
1280 4244 : do irhoij=1,nrhoij_in
1281 2062 : cplex =pawrhoij_in(irhoij)%cplex_rhoij
1282 2062 : qphase =pawrhoij_in(irhoij)%qphase
1283 2062 : lmn2_size=pawrhoij_in(irhoij)%lmn2_size
1284 2062 : nspden =pawrhoij_in(irhoij)%nspden
1285 2062 : if (with_lmnmix_) lmnmix=pawrhoij_in(irhoij)%lmnmix_sz
1286 2062 : if (with_grhoij_) ngrhoij=pawrhoij_in(irhoij)%ngrhoij
1287 2062 : if (with_rhoijp_) use_rhoijp=pawrhoij_in(irhoij)%use_rhoijp
1288 2062 : if (with_rhoijres_)use_rhoijres=pawrhoij_in(irhoij)%use_rhoijres
1289 2062 : if (with_rhoij__) use_rhoij_=pawrhoij_in(irhoij)%use_rhoij_
1290 2062 : buf_int_size=buf_int_size+16
1291 2062 : if (use_rhoijp>0) then
1292 2062 : nselect=pawrhoij_in(irhoij)%nrhoijsel
1293 2062 : buf_int_size=buf_int_size+nselect
1294 2062 : buf_dp_size=buf_dp_size + cplex*qphase*nselect*nspden
1295 : end if
1296 2062 : if (lmnmix>0) buf_int_size=buf_int_size+lmnmix
1297 2062 : if (ngrhoij>0) buf_dp_size=buf_dp_size + cplex*qphase*lmn2_size*nspden*ngrhoij
1298 2062 : if (use_rhoijres>0) buf_dp_size=buf_dp_size + cplex*qphase*lmn2_size*nspden
1299 4244 : if (use_rhoij_>0) then
1300 0 : rhoij_size2=size(pawrhoij_in(irhoij)%rhoij_,dim=2)
1301 0 : buf_dp_size=buf_dp_size + cplex*qphase*lmn2_size*rhoij_size2
1302 : end if
1303 : end do
1304 :
1305 : !Fill input buffers
1306 6546 : LIBPAW_ALLOCATE(buf_int,(buf_int_size))
1307 6546 : LIBPAW_ALLOCATE(buf_dp ,(buf_dp_size))
1308 2182 : indx_int=1;indx_dp =1
1309 2182 : lmnmix=0;ngrhoij=0;nselect=0;rhoij_size2=0
1310 2182 : use_rhoijp=0;use_rhoijres=0;use_rhoij_=0
1311 4244 : do irhoij=1,nrhoij_in
1312 2062 : cplex =pawrhoij_in(irhoij)%cplex_rhoij
1313 2062 : qphase =pawrhoij_in(irhoij)%qphase
1314 2062 : lmn2_size=pawrhoij_in(irhoij)%lmn2_size
1315 2062 : nspden =pawrhoij_in(irhoij)%nspden
1316 2062 : if (with_lmnmix_) lmnmix=pawrhoij_in(irhoij)%lmnmix_sz
1317 2062 : if (with_grhoij_) ngrhoij=pawrhoij_in(irhoij)%ngrhoij
1318 2062 : if (with_rhoijp_) use_rhoijp=pawrhoij_in(irhoij)%use_rhoijp
1319 2062 : if (use_rhoijp > 0) nselect=pawrhoij_in(irhoij)%nrhoijsel
1320 2062 : if (with_rhoijres_)use_rhoijres=pawrhoij_in(irhoij)%use_rhoijres
1321 2062 : if (with_rhoij__) use_rhoij_ =pawrhoij_in(irhoij)%use_rhoij_
1322 2062 : if (use_rhoij_> 0) rhoij_size2 =size(pawrhoij_in(irhoij)%rhoij_,dim=2)
1323 2062 : buf_int(indx_int)=my_atmtab(irhoij) ;indx_int=indx_int+1
1324 2062 : buf_int(indx_int)=cplex ;indx_int=indx_int+1
1325 2062 : buf_int(indx_int)=qphase ;indx_int=indx_int+1
1326 2062 : buf_int(indx_int)=lmn2_size ;indx_int=indx_int+1
1327 2062 : buf_int(indx_int)=nspden ;indx_int=indx_int+1
1328 2062 : buf_int(indx_int)=nselect ;indx_int=indx_int+1
1329 2062 : buf_int(indx_int)=lmnmix ;indx_int=indx_int+1
1330 2062 : buf_int(indx_int)=ngrhoij ;indx_int=indx_int+1
1331 2062 : buf_int(indx_int)=use_rhoijp ;indx_int=indx_int+1
1332 2062 : buf_int(indx_int)=use_rhoijres ;indx_int=indx_int+1
1333 2062 : buf_int(indx_int)=use_rhoij_ ;indx_int=indx_int+1
1334 2062 : buf_int(indx_int)=rhoij_size2 ;indx_int=indx_int+1
1335 2062 : buf_int(indx_int)=pawrhoij_in(irhoij)%itypat ;indx_int=indx_int+1
1336 2062 : buf_int(indx_int)=pawrhoij_in(irhoij)%lmn_size ;indx_int=indx_int+1
1337 2062 : buf_int(indx_int)=pawrhoij_in(irhoij)%nsppol ;indx_int=indx_int+1
1338 2062 : buf_int(indx_int)=pawrhoij_in(irhoij)%nspinor ;indx_int=indx_int+1
1339 2062 : if (use_rhoijp>0) then
1340 82026 : buf_int(indx_int:indx_int+nselect-1)=pawrhoij_in(irhoij)%rhoijselect(1:nselect)
1341 4620 : indx_int=indx_int+nselect
1342 4620 : do isp=1,nspden
1343 7214 : do ii=1,qphase
1344 2594 : jj=(ii-1)*cplex*lmn2_size
1345 : buf_dp(indx_dp:indx_dp+cplex*nselect-1)= &
1346 140032 : & pawrhoij_in(irhoij)%rhoijp(jj+1:jj+cplex*nselect,isp)
1347 5152 : indx_dp=indx_dp+cplex*nselect
1348 : end do
1349 : end do
1350 : end if
1351 2062 : if (lmnmix>0) then
1352 36758 : buf_int(indx_int:indx_int+lmnmix-1)=pawrhoij_in(irhoij)%kpawmix(1:lmnmix)
1353 : indx_int=indx_int+lmnmix
1354 : end if
1355 2062 : if (ngrhoij>0) then
1356 0 : do isp=1,nspden
1357 0 : do ii=1,cplex*qphase*lmn2_size
1358 0 : buf_dp(indx_dp:indx_dp+ngrhoij-1)=pawrhoij_in(irhoij)%grhoij(1:ngrhoij,ii,isp)
1359 0 : indx_dp=indx_dp+ngrhoij
1360 : end do
1361 : end do
1362 : end if
1363 2062 : if (use_rhoijres>0) then
1364 1146 : do isp=1,nspden
1365 : buf_dp(indx_dp:indx_dp+cplex*qphase*lmn2_size-1)= &
1366 58404 : pawrhoij_in(irhoij)%rhoijres(1:cplex*qphase*lmn2_size,isp)
1367 1146 : indx_dp=indx_dp+cplex*qphase*lmn2_size
1368 : end do
1369 : end if
1370 4244 : if (use_rhoij_>0) then
1371 0 : do isp=1,rhoij_size2
1372 : buf_dp(indx_dp:indx_dp+cplex*qphase*lmn2_size-1)= &
1373 0 : & pawrhoij_in(irhoij)%rhoij_(1:cplex*qphase*lmn2_size,isp)
1374 0 : indx_dp=indx_dp+cplex*qphase*lmn2_size
1375 : end do
1376 : end if
1377 : end do
1378 :
1379 : !Check
1380 2182 : if ((indx_int-1/=buf_int_size).or.(indx_dp-1/=buf_dp_size)) then
1381 0 : write(msg,*) 'Wrong buffer sizes: buf_int_size=',buf_int_size,' buf_dp_size=',buf_dp_size
1382 0 : LIBPAW_BUG(msg)
1383 : end if
1384 :
1385 : !Communicate (1 gather for integers, 1 gather for reals)
1386 6546 : LIBPAW_ALLOCATE(count_int,(nproc_atom))
1387 4364 : LIBPAW_ALLOCATE(displ_int,(nproc_atom))
1388 4364 : LIBPAW_ALLOCATE(count_dp ,(nproc_atom))
1389 4364 : LIBPAW_ALLOCATE(displ_dp ,(nproc_atom))
1390 6546 : LIBPAW_ALLOCATE(count_tot,(2*nproc_atom))
1391 2182 : bufsz(1)=buf_int_size; bufsz(2)=buf_dp_size
1392 2182 : call xmpi_allgather(bufsz,2,count_tot,comm_atom,ierr)
1393 9554 : do ii=1,nproc_atom
1394 7372 : count_int(ii)=count_tot(2*ii-1)
1395 9554 : count_dp (ii)=count_tot(2*ii)
1396 : end do
1397 2182 : displ_int(1)=0;displ_dp(1)=0
1398 7372 : do ii=2,nproc_atom
1399 5190 : displ_int(ii)=displ_int(ii-1)+count_int(ii-1)
1400 7372 : displ_dp (ii)=displ_dp (ii-1)+count_dp (ii-1)
1401 : end do
1402 9554 : buf_int_size_all=sum(count_int)
1403 9554 : buf_dp_size_all =sum(count_dp)
1404 2182 : LIBPAW_DEALLOCATE(count_tot)
1405 2182 : if (master==-1.or.me_atom==master) then
1406 5637 : LIBPAW_ALLOCATE(buf_int_all,(buf_int_size_all))
1407 5637 : LIBPAW_ALLOCATE(buf_dp_all ,(buf_dp_size_all))
1408 : else
1409 303 : LIBPAW_ALLOCATE(buf_int_all,(0))
1410 303 : LIBPAW_ALLOCATE(buf_dp_all ,(0))
1411 : end if
1412 2182 : if (master==-1) then
1413 1740 : call xmpi_allgatherv(buf_int,buf_int_size,buf_int_all,count_int,displ_int,comm_atom,ierr)
1414 1740 : call xmpi_allgatherv(buf_dp ,buf_dp_size ,buf_dp_all ,count_dp ,displ_dp ,comm_atom,ierr)
1415 : else
1416 442 : call xmpi_gatherv(buf_int,buf_int_size,buf_int_all,count_int,displ_int,master,comm_atom,ierr)
1417 442 : call xmpi_gatherv(buf_dp ,buf_dp_size ,buf_dp_all ,count_dp ,displ_dp ,master,comm_atom,ierr)
1418 : end if
1419 2182 : LIBPAW_DEALLOCATE(count_int)
1420 2182 : LIBPAW_DEALLOCATE(displ_int)
1421 2182 : LIBPAW_DEALLOCATE(count_dp)
1422 2182 : LIBPAW_DEALLOCATE(displ_dp)
1423 :
1424 : !Retrieve data from output buffer
1425 2182 : if (master==-1.or.me_atom==master) then
1426 1879 : indx_int=1;indx_dp=1
1427 1879 : call pawrhoij_free(pawrhoij_gathered)
1428 7523 : do irhoij=1,nrhoij_out
1429 5644 : jrhoij =buf_int_all(indx_int) ;indx_int=indx_int+1
1430 5644 : cplex =buf_int_all(indx_int) ;indx_int=indx_int+1
1431 5644 : qphase =buf_int_all(indx_int) ;indx_int=indx_int+1
1432 5644 : lmn2_size =buf_int_all(indx_int) ;indx_int=indx_int+1
1433 5644 : nspden =buf_int_all(indx_int) ;indx_int=indx_int+1
1434 5644 : nselect =buf_int_all(indx_int) ;indx_int=indx_int+1
1435 5644 : lmnmix =buf_int_all(indx_int) ;indx_int=indx_int+1
1436 5644 : ngrhoij =buf_int_all(indx_int) ;indx_int=indx_int+1
1437 5644 : use_rhoijp =buf_int_all(indx_int) ;indx_int=indx_int+1
1438 5644 : use_rhoijres=buf_int_all(indx_int) ;indx_int=indx_int+1
1439 5644 : use_rhoij_ =buf_int_all(indx_int) ;indx_int=indx_int+1
1440 5644 : rhoij_size2 =buf_int_all(indx_int) ;indx_int=indx_int+1
1441 5644 : pawrhoij_gathered(jrhoij)%itypat=buf_int_all(indx_int) ;indx_int=indx_int+1
1442 5644 : pawrhoij_gathered(jrhoij)%lmn_size=buf_int_all(indx_int) ;indx_int=indx_int+1
1443 5644 : pawrhoij_gathered(jrhoij)%nsppol=buf_int_all(indx_int) ;indx_int=indx_int+1
1444 5644 : pawrhoij_gathered(jrhoij)%nspinor=buf_int_all(indx_int) ;indx_int=indx_int+1
1445 5644 : pawrhoij_gathered(jrhoij)%cplex_rhoij=cplex
1446 5644 : pawrhoij_gathered(jrhoij)%qphase=qphase
1447 5644 : pawrhoij_gathered(jrhoij)%lmn2_size=lmn2_size
1448 5644 : pawrhoij_gathered(jrhoij)%nspden=nspden
1449 5644 : pawrhoij_gathered(jrhoij)%nrhoijsel=nselect
1450 5644 : pawrhoij_gathered(jrhoij)%lmnmix_sz=lmnmix
1451 5644 : pawrhoij_gathered(jrhoij)%ngrhoij=ngrhoij
1452 5644 : pawrhoij_gathered(jrhoij)%use_rhoijp=use_rhoijp
1453 5644 : pawrhoij_gathered(jrhoij)%use_rhoijres=use_rhoijres
1454 5644 : pawrhoij_gathered(jrhoij)%use_rhoij_=use_rhoij_
1455 5644 : if (use_rhoijp>0) then
1456 16932 : LIBPAW_ALLOCATE(pawrhoij_gathered(jrhoij)%rhoijselect,(lmn2_size))
1457 208234 : pawrhoij_gathered(jrhoij)%rhoijselect(1:nselect)=buf_int_all(indx_int:indx_int+nselect-1)
1458 192494 : if (nselect < lmn2_size )pawrhoij_gathered(jrhoij)%rhoijselect(nselect+1:lmn2_size)=0
1459 5644 : indx_int=indx_int+nselect
1460 22576 : LIBPAW_ALLOCATE(pawrhoij_gathered(jrhoij)%rhoijp,(qphase*cplex*lmn2_size,nspden))
1461 12666 : do isp=1,nspden
1462 19768 : do ii=1,qphase
1463 7102 : jj=(ii-1)*cplex*lmn2_size
1464 : pawrhoij_gathered(jrhoij)%rhoijp(jj+1:jj+cplex*nselect,isp)= &
1465 337506 : & buf_dp_all(indx_dp:indx_dp+cplex*nselect-1)
1466 7102 : if (nselect<lmn2_size) &
1467 249588 : & pawrhoij_gathered(jrhoij)%rhoijp(jj+cplex*nselect+1:jj+cplex*lmn2_size,isp)=zero
1468 14124 : indx_dp=indx_dp+cplex*nselect
1469 : end do
1470 : end do
1471 : end if
1472 5644 : if (lmnmix>0) then
1473 4560 : LIBPAW_ALLOCATE(pawrhoij_gathered(jrhoij)%kpawmix,(lmnmix))
1474 116208 : pawrhoij_gathered(jrhoij)%kpawmix(1:lmnmix)=buf_int_all(indx_int:indx_int+lmnmix-1)
1475 : indx_int=indx_int+lmnmix
1476 : end if
1477 5644 : if (ngrhoij>0) then
1478 0 : LIBPAW_ALLOCATE(pawrhoij_gathered(jrhoij)%grhoij,(ngrhoij,qphase*cplex*lmn2_size,nspden))
1479 0 : do isp=1,nspden
1480 0 : do ii=1,cplex*qphase*lmn2_size
1481 0 : pawrhoij_gathered(jrhoij)%grhoij(1:ngrhoij,ii,isp)=buf_dp_all(indx_dp:indx_dp+ngrhoij-1)
1482 0 : indx_dp=indx_dp+ngrhoij
1483 : end do
1484 : end do
1485 : end if
1486 5644 : if (use_rhoijres>0) then
1487 6080 : LIBPAW_ALLOCATE(pawrhoij_gathered(jrhoij)%rhoijres,(qphase*cplex*lmn2_size,nspden))
1488 3464 : do isp=1,nspden
1489 : pawrhoij_gathered(jrhoij)%rhoijres(1:cplex*qphase*lmn2_size,isp)= &
1490 173104 : & buf_dp_all(indx_dp:indx_dp+cplex*qphase*lmn2_size-1)
1491 3464 : indx_dp=indx_dp+cplex*qphase*lmn2_size
1492 : end do
1493 : end if
1494 7523 : if (use_rhoij_>0) then
1495 0 : LIBPAW_ALLOCATE(pawrhoij_gathered(jrhoij)%rhoij_,(qphase*cplex*lmn2_size,rhoij_size2))
1496 0 : do isp=1,rhoij_size2
1497 : pawrhoij_gathered(jrhoij)%rhoij_(1:cplex*qphase*lmn2_size,isp)= &
1498 0 : & buf_dp_all(indx_dp:indx_dp+cplex*qphase*lmn2_size-1)
1499 0 : indx_dp=indx_dp+cplex*qphase*lmn2_size
1500 : end do
1501 : end if
1502 : end do
1503 1879 : if ((indx_int/=1+buf_int_size_all).or.(indx_dp/=1+buf_dp_size_all)) then
1504 0 : write(msg,*) 'Wrong buffer sizes: buf_int_size_all=',buf_int_size_all,' buf_dp_size_all=',buf_dp_size_all
1505 0 : LIBPAW_BUG(msg)
1506 : end if
1507 : end if
1508 :
1509 : !Free memory
1510 2182 : call free_my_atmtab(my_atmtab,my_atmtab_allocated)
1511 2182 : LIBPAW_DEALLOCATE(buf_int)
1512 2182 : LIBPAW_DEALLOCATE(buf_dp)
1513 2182 : LIBPAW_DEALLOCATE(buf_int_all)
1514 2182 : LIBPAW_DEALLOCATE(buf_dp_all)
1515 :
1516 6546 : end subroutine pawrhoij_gather
1517 : !!***
1518 :
1519 : !----------------------------------------------------------------------
1520 :
1521 : !!****f* m_pawrhoij/pawrhoij_bcast
1522 : !! NAME
1523 : !! pawrhoij_bcast
1524 : !!
1525 : !! FUNCTION
1526 : !! Broadcast pawrhoij datastructures
1527 : !! Can take into account a distribution of data over a "atom" communicator
1528 : !!
1529 : !! INPUTS
1530 : !! master=master communicator receiving data
1531 : !! mpicomm= MPI communicator
1532 : !! comm_atom= --optional-- MPI communicator over atoms
1533 : !! pawrhoij_in(:)<type(pawrhoij_type)>= input rhoij datastructures on master process
1534 : !!
1535 : !! OUTPUT
1536 : !! pawrhoij_out(:)<type(pawrhoij_type)>= output rhoij datastructure on every process
1537 : !! Eventually distributed according to comm_atom communicator
1538 : !!
1539 : !! SOURCE
1540 :
1541 0 : subroutine pawrhoij_bcast(pawrhoij_in,pawrhoij_out,master,mpicomm,comm_atom)
1542 :
1543 : !Arguments ------------------------------------
1544 : !scalars
1545 : integer,intent(in) :: master,mpicomm
1546 : integer,intent(in),optional :: comm_atom
1547 : !arrays
1548 : type(pawrhoij_type),intent(in) :: pawrhoij_in(:)
1549 : type(pawrhoij_type),intent(inout) :: pawrhoij_out(:)
1550 : !Local variables-------------------------------
1551 : !scalars
1552 : integer :: buf_dp_size,buf_dp_size_all,buf_int_size,buf_int_size_all
1553 : integer :: cplex,ierr,ii,indx_dp,indx_int,iproc,irhoij,isp,jj,jrhoij,lmn2_size,lmnmix,me,me_atom
1554 : integer :: my_comm_atom,ngrhoij,nproc,nproc_atom,nrhoij_in,nrhoij_out,nrhoij_out_all
1555 : integer :: nselect,nspden,qphase,rhoij_size2,use_rhoijp,use_rhoijres,use_rhoij_
1556 : logical :: my_atmtab_allocated,paral_atom
1557 : character(len=500) :: msg
1558 : !arrays
1559 : integer :: buf_size(2)
1560 0 : integer,allocatable :: atmtab(:),buf_int_size_i(:),buf_dp_size_i(:)
1561 0 : integer,allocatable :: count_dp(:),count_int(:),disp_dp(:),disp_int(:),nrhoij_out_i(:)
1562 0 : integer,allocatable,target :: buf_int(:)
1563 0 : integer,pointer :: buf_int_all(:),my_atmtab(:)
1564 0 : real(dp),allocatable,target :: buf_dp(:)
1565 0 : real(dp),pointer :: buf_dp_all(:)
1566 :
1567 : ! *************************************************************************
1568 :
1569 : !Load MPI "atom" distribution data
1570 0 : my_comm_atom=xmpi_comm_self;nproc_atom=1;me_atom=0
1571 0 : if (present(comm_atom)) then
1572 0 : my_comm_atom=comm_atom
1573 0 : me_atom=xmpi_comm_rank(my_comm_atom)
1574 0 : nproc_atom=xmpi_comm_size(my_comm_atom)
1575 0 : paral_atom=(nproc_atom>1)
1576 0 : if (my_comm_atom/=mpicomm.and.nproc_atom/=1) then
1577 0 : msg='wrong comm_atom communicator !'
1578 0 : LIBPAW_BUG(msg)
1579 : end if
1580 : end if
1581 :
1582 : !Load global MPI data
1583 0 : me=xmpi_comm_rank(mpicomm)
1584 0 : nproc=xmpi_comm_size(mpicomm)
1585 :
1586 : !Just copy in case of a sequential run
1587 0 : if (nproc==1.and.nproc_atom==1) then
1588 0 : call pawrhoij_copy(pawrhoij_in,pawrhoij_out,.false.,.false.,.false.)
1589 0 : return
1590 : end if
1591 :
1592 : !Retrieve and test pawrhoij sizes
1593 0 : nrhoij_in=0;if (me==master) nrhoij_in=size(pawrhoij_in)
1594 0 : nrhoij_out=size(pawrhoij_out);nrhoij_out_all=nrhoij_out
1595 0 : if (paral_atom) then
1596 0 : LIBPAW_ALLOCATE(nrhoij_out_i,(nproc_atom))
1597 0 : buf_size(1)=nrhoij_out
1598 0 : call xmpi_allgather(buf_size,1,nrhoij_out_i,my_comm_atom,ierr)
1599 0 : nrhoij_out_all=sum(nrhoij_out_i)
1600 : end if
1601 0 : if (me==master.and.nrhoij_in/=nrhoij_out_all) then
1602 0 : msg='pawrhoij_in or pawrhoij_out wrongly allocated!'
1603 0 : LIBPAW_BUG(msg)
1604 : end if
1605 :
1606 : !Retrieve table(s) of atoms (if necessary)
1607 0 : if (paral_atom) then
1608 0 : nullify(my_atmtab)
1609 : call get_my_atmtab(my_comm_atom,my_atmtab,my_atmtab_allocated,paral_atom, &
1610 0 : & nrhoij_out_all,my_natom_ref=nrhoij_out)
1611 0 : LIBPAW_ALLOCATE(disp_int,(nproc_atom))
1612 0 : disp_int(1)=0
1613 0 : do iproc=2,nproc_atom
1614 0 : disp_int(iproc)=disp_int(iproc-1)+nrhoij_out_i(iproc-1)
1615 : end do
1616 0 : LIBPAW_ALLOCATE(atmtab,(nrhoij_in))
1617 : call xmpi_gatherv(my_atmtab,nrhoij_out,atmtab,nrhoij_out_i,disp_int,&
1618 0 : & master,my_comm_atom,ierr)
1619 0 : LIBPAW_DEALLOCATE(disp_int)
1620 : end if
1621 :
1622 : !Compute sizes of input buffers and broadcast them
1623 0 : LIBPAW_ALLOCATE(buf_int_size_i,(nrhoij_out_all))
1624 0 : LIBPAW_ALLOCATE(buf_dp_size_i ,(nrhoij_out_all))
1625 0 : if (me==master) then
1626 0 : buf_int_size_i(:)=0;buf_dp_size_i(:)=0
1627 0 : do irhoij=1,nrhoij_in
1628 0 : jrhoij=irhoij;if (paral_atom) jrhoij=atmtab(irhoij)
1629 0 : cplex =pawrhoij_in(jrhoij)%cplex_rhoij
1630 0 : qphase =pawrhoij_in(jrhoij)%qphase
1631 0 : lmn2_size =pawrhoij_in(jrhoij)%lmn2_size
1632 0 : nspden =pawrhoij_in(jrhoij)%nspden
1633 0 : lmnmix =pawrhoij_in(jrhoij)%lmnmix_sz
1634 0 : ngrhoij =pawrhoij_in(jrhoij)%ngrhoij
1635 0 : use_rhoijp =pawrhoij_in(jrhoij)%use_rhoijp
1636 0 : use_rhoijres=pawrhoij_in(jrhoij)%use_rhoijres
1637 0 : use_rhoij_ =pawrhoij_in(jrhoij)%use_rhoij_
1638 0 : buf_int_size_i(irhoij)=buf_int_size_i(irhoij)+16
1639 0 : if (ngrhoij>0) buf_dp_size_i(irhoij)=buf_dp_size_i(irhoij)+cplex*qphase*lmn2_size*nspden*ngrhoij
1640 0 : if (use_rhoijres>0) buf_dp_size_i(irhoij)=buf_dp_size_i(irhoij)+cplex*qphase*lmn2_size*nspden
1641 0 : if (use_rhoijp>0) then
1642 0 : nselect=pawrhoij_in(jrhoij)%nrhoijsel
1643 0 : buf_int_size_i(irhoij)=buf_int_size_i(irhoij)+nselect
1644 0 : buf_dp_size_i(irhoij)=buf_dp_size_i(irhoij)+cplex*qphase*nselect*nspden
1645 : end if
1646 0 : if (use_rhoij_>0) then
1647 0 : rhoij_size2=size(pawrhoij_in(jrhoij)%rhoij_,dim=2)
1648 0 : buf_dp_size_i(irhoij)=buf_dp_size_i(irhoij)+cplex*qphase*lmn2_size*rhoij_size2
1649 : end if
1650 : end do
1651 : end if
1652 0 : call xmpi_bcast(buf_int_size_i,master,mpicomm,ierr)
1653 0 : call xmpi_bcast(buf_dp_size_i,master,mpicomm,ierr)
1654 0 : buf_int_size_all=sum(buf_int_size_i) ; buf_dp_size_all=sum(buf_dp_size_i)
1655 :
1656 : !Prepare buffers/tabs for communication
1657 0 : if (paral_atom) then
1658 0 : LIBPAW_ALLOCATE(count_int,(nproc_atom))
1659 0 : LIBPAW_ALLOCATE(count_dp,(nproc_atom))
1660 0 : LIBPAW_ALLOCATE(disp_int,(nproc_atom))
1661 0 : LIBPAW_ALLOCATE(disp_dp,(nproc_atom))
1662 0 : indx_int=0
1663 0 : do iproc=1,nproc_atom
1664 0 : ii=nrhoij_out_i(iproc)
1665 0 : count_int(iproc)=sum(buf_int_size_i(indx_int+1:indx_int+ii))
1666 0 : count_dp (iproc)=sum(buf_dp_size_i (indx_int+1:indx_int+ii))
1667 0 : indx_int=indx_int+ii
1668 : end do
1669 0 : disp_int(1)=0;disp_dp(1)=0
1670 0 : do iproc=2,nproc_atom
1671 0 : disp_int(iproc)=disp_int(iproc-1)+count_int(iproc-1)
1672 0 : disp_dp (iproc)=disp_dp (iproc-1)+count_dp (iproc-1)
1673 : end do
1674 0 : if (buf_int_size_all/=sum(count_int).or.buf_dp_size_all/=sum(count_dp)) then
1675 0 : msg='(1) Wrong buffer sizes !'
1676 0 : LIBPAW_BUG(msg)
1677 : end if
1678 0 : buf_int_size=count_int(me_atom+1)
1679 0 : buf_dp_size =count_dp(me_atom+1)
1680 0 : LIBPAW_ALLOCATE(buf_int,(buf_int_size))
1681 0 : LIBPAW_ALLOCATE(buf_dp ,(buf_dp_size))
1682 0 : if (me==master) then
1683 0 : LIBPAW_POINTER_ALLOCATE(buf_int_all,(buf_int_size_all))
1684 0 : LIBPAW_POINTER_ALLOCATE(buf_dp_all ,(buf_dp_size_all))
1685 : else
1686 0 : LIBPAW_POINTER_ALLOCATE(buf_int_all,(1))
1687 0 : LIBPAW_POINTER_ALLOCATE(buf_dp_all ,(1))
1688 : end if
1689 0 : LIBPAW_DEALLOCATE(nrhoij_out_i)
1690 : else
1691 0 : buf_int_size=buf_int_size_all
1692 0 : buf_dp_size =buf_dp_size_all
1693 0 : LIBPAW_ALLOCATE(buf_int,(buf_int_size))
1694 0 : LIBPAW_ALLOCATE(buf_dp ,(buf_dp_size))
1695 0 : buf_int_all => buf_int
1696 0 : buf_dp_all => buf_dp
1697 : end if
1698 0 : LIBPAW_DEALLOCATE(buf_int_size_i)
1699 0 : LIBPAW_DEALLOCATE(buf_dp_size_i)
1700 :
1701 : !Fill input buffers
1702 0 : if (me==master) then
1703 : indx_int=1;indx_dp =1
1704 0 : do irhoij=1,nrhoij_in
1705 0 : jrhoij=irhoij;if (paral_atom) jrhoij=atmtab(irhoij)
1706 0 : cplex =pawrhoij_in(jrhoij)%cplex_rhoij
1707 0 : qphase =pawrhoij_in(jrhoij)%qphase
1708 0 : lmn2_size =pawrhoij_in(jrhoij)%lmn2_size
1709 0 : nspden =pawrhoij_in(jrhoij)%nspden
1710 0 : lmnmix =pawrhoij_in(jrhoij)%lmnmix_sz
1711 0 : ngrhoij =pawrhoij_in(jrhoij)%ngrhoij
1712 0 : use_rhoijp =pawrhoij_in(jrhoij)%use_rhoijp
1713 0 : nselect =pawrhoij_in(jrhoij)%nrhoijsel
1714 0 : use_rhoijres=pawrhoij_in(jrhoij)%use_rhoijres
1715 0 : use_rhoij_ =pawrhoij_in(jrhoij)%use_rhoij_
1716 0 : rhoij_size2 =size(pawrhoij_in(jrhoij)%rhoij_,dim=2)
1717 0 : buf_int_all(indx_int)=jrhoij ;indx_int=indx_int+1 ! Not used !
1718 0 : buf_int_all(indx_int)=cplex ;indx_int=indx_int+1
1719 0 : buf_int_all(indx_int)=qphase ;indx_int=indx_int+1
1720 0 : buf_int_all(indx_int)=lmn2_size ;indx_int=indx_int+1
1721 0 : buf_int_all(indx_int)=nspden ;indx_int=indx_int+1
1722 0 : buf_int_all(indx_int)=nselect ;indx_int=indx_int+1
1723 0 : buf_int_all(indx_int)=lmnmix ;indx_int=indx_int+1
1724 0 : buf_int_all(indx_int)=ngrhoij ;indx_int=indx_int+1
1725 0 : buf_int_all(indx_int)=use_rhoijp ;indx_int=indx_int+1
1726 0 : buf_int_all(indx_int)=use_rhoijres ;indx_int=indx_int+1
1727 0 : buf_int_all(indx_int)=use_rhoij_ ;indx_int=indx_int+1
1728 0 : buf_int_all(indx_int)=rhoij_size2 ;indx_int=indx_int+1
1729 0 : buf_int_all(indx_int)=pawrhoij_in(jrhoij)%itypat ;indx_int=indx_int+1
1730 0 : buf_int_all(indx_int)=pawrhoij_in(jrhoij)%lmn_size;indx_int=indx_int+1
1731 0 : buf_int_all(indx_int)=pawrhoij_in(jrhoij)%nsppol ;indx_int=indx_int+1
1732 0 : buf_int_all(indx_int)=pawrhoij_in(jrhoij)%nspinor ;indx_int=indx_int+1
1733 0 : if (use_rhoijp>0) then
1734 0 : buf_int_all(indx_int:indx_int+nselect-1)=pawrhoij_in(jrhoij)%rhoijselect(1:nselect)
1735 0 : indx_int=indx_int+nselect
1736 0 : do isp=1,nspden
1737 0 : do ii=1,qphase
1738 0 : jj=(ii-1)*cplex*lmn2_size
1739 : buf_dp_all(indx_dp:indx_dp+cplex*nselect-1)= &
1740 0 : & pawrhoij_in(jrhoij)%rhoijp(jj+1:jj+cplex*nselect,isp)
1741 0 : indx_dp=indx_dp+cplex*nselect
1742 : end do
1743 : end do
1744 : end if
1745 0 : if (lmnmix>0) then
1746 0 : buf_int_all(indx_int:indx_int+lmnmix-1)=pawrhoij_in(jrhoij)%kpawmix(1:lmnmix)
1747 : indx_int=indx_int+lmnmix
1748 : end if
1749 0 : if (ngrhoij>0) then
1750 0 : do isp=1,nspden
1751 0 : do ii=1,cplex*qphase*lmn2_size
1752 0 : buf_dp_all(indx_dp:indx_dp+ngrhoij-1)=pawrhoij_in(jrhoij)%grhoij(1:ngrhoij,ii,isp)
1753 0 : indx_dp=indx_dp+ngrhoij
1754 : end do
1755 : end do
1756 : end if
1757 0 : if (use_rhoijres>0) then
1758 0 : do isp=1,nspden
1759 : buf_dp_all(indx_dp:indx_dp+cplex*qphase*lmn2_size-1)= &
1760 0 : & pawrhoij_in(jrhoij)%rhoijres(1:cplex*qphase*lmn2_size,isp)
1761 0 : indx_dp=indx_dp+cplex*qphase*lmn2_size
1762 : end do
1763 : end if
1764 0 : if (use_rhoij_>0) then
1765 0 : do isp=1,rhoij_size2
1766 : buf_dp_all(indx_dp:indx_dp+cplex*qphase*lmn2_size-1)= &
1767 0 : & pawrhoij_in(jrhoij)%rhoij_(1:cplex*qphase*lmn2_size,isp)
1768 0 : indx_dp=indx_dp+cplex*qphase*lmn2_size
1769 : end do
1770 : end if
1771 : end do
1772 : ! Check
1773 0 : if ((indx_int-1/=buf_int_size_all).or.(indx_dp-1/=buf_dp_size_all)) then
1774 0 : msg='(2) Wrong buffer sizes !'
1775 0 : LIBPAW_BUG(msg)
1776 : end if
1777 : end if ! me=master
1778 :
1779 : !Communicate
1780 0 : if (paral_atom) then
1781 0 : call xmpi_scatterv(buf_int_all,count_int,disp_int,buf_int,buf_int_size,master,mpicomm,ierr)
1782 0 : call xmpi_scatterv(buf_dp_all ,count_dp ,disp_dp ,buf_dp ,buf_dp_size ,master,mpicomm,ierr)
1783 : else
1784 0 : call xmpi_bcast(buf_int,master,mpicomm,ierr)
1785 0 : call xmpi_bcast(buf_dp ,master,mpicomm,ierr)
1786 : end if
1787 :
1788 : !Retrieve data from output buffer
1789 0 : indx_int=1;indx_dp=1
1790 0 : call pawrhoij_free(pawrhoij_out)
1791 0 : do irhoij=1,nrhoij_out
1792 0 : jrhoij =buf_int(indx_int);indx_int=indx_int+1 ! Not used !
1793 0 : cplex =buf_int(indx_int);indx_int=indx_int+1
1794 0 : qphase =buf_int(indx_int);indx_int=indx_int+1
1795 0 : lmn2_size =buf_int(indx_int);indx_int=indx_int+1
1796 0 : nspden =buf_int(indx_int);indx_int=indx_int+1
1797 0 : nselect =buf_int(indx_int);indx_int=indx_int+1
1798 0 : lmnmix =buf_int(indx_int);indx_int=indx_int+1
1799 0 : ngrhoij =buf_int(indx_int);indx_int=indx_int+1
1800 0 : use_rhoijp =buf_int(indx_int);indx_int=indx_int+1
1801 0 : use_rhoijres=buf_int(indx_int);indx_int=indx_int+1
1802 0 : use_rhoij_ =buf_int(indx_int);indx_int=indx_int+1
1803 0 : rhoij_size2 =buf_int(indx_int);indx_int=indx_int+1
1804 0 : pawrhoij_out(irhoij)%itypat=buf_int(indx_int) ;indx_int=indx_int+1
1805 0 : pawrhoij_out(irhoij)%lmn_size=buf_int(indx_int);indx_int=indx_int+1
1806 0 : pawrhoij_out(irhoij)%nsppol=buf_int(indx_int) ;indx_int=indx_int+1
1807 0 : pawrhoij_out(irhoij)%nspinor=buf_int(indx_int) ;indx_int=indx_int+1
1808 0 : pawrhoij_out(irhoij)%cplex_rhoij=cplex
1809 0 : pawrhoij_out(irhoij)%qphase=qphase
1810 0 : pawrhoij_out(irhoij)%lmn2_size=lmn2_size
1811 0 : pawrhoij_out(irhoij)%nspden=nspden
1812 0 : pawrhoij_out(irhoij)%nrhoijsel=nselect
1813 0 : pawrhoij_out(irhoij)%lmnmix_sz=lmnmix
1814 0 : pawrhoij_out(irhoij)%ngrhoij=ngrhoij
1815 0 : pawrhoij_out(irhoij)%use_rhoijp=use_rhoijp
1816 0 : pawrhoij_out(irhoij)%use_rhoijres=use_rhoijres
1817 0 : pawrhoij_out(irhoij)%use_rhoij_=use_rhoij_
1818 0 : if (use_rhoijp>0) then
1819 0 : LIBPAW_ALLOCATE(pawrhoij_out(irhoij)%rhoijselect,(nselect))
1820 0 : pawrhoij_out(irhoij)%rhoijselect=0
1821 0 : pawrhoij_out(irhoij)%rhoijselect(1:nselect)=buf_int(indx_int:indx_int+nselect-1)
1822 0 : indx_int=indx_int+nselect
1823 0 : LIBPAW_ALLOCATE(pawrhoij_out(irhoij)%rhoijp,(qphase*cplex*nselect,nspden))
1824 0 : do isp=1,nspden
1825 0 : do ii=1,qphase
1826 0 : jj=(ii-1)*cplex*lmn2_size
1827 : pawrhoij_out(irhoij)%rhoijp(jj+1:jj+cplex*nselect,isp)= &
1828 0 : & buf_dp(indx_dp:indx_dp+cplex*nselect-1)
1829 0 : indx_dp=indx_dp+cplex*nselect
1830 : end do
1831 : end do
1832 : end if
1833 0 : if (lmnmix>0) then
1834 0 : LIBPAW_ALLOCATE(pawrhoij_out(irhoij)%kpawmix,(lmnmix))
1835 0 : pawrhoij_out(irhoij)%kpawmix(1:lmnmix)=buf_int(indx_int:indx_int+lmnmix-1)
1836 : indx_int=indx_int+lmnmix
1837 : end if
1838 0 : if (ngrhoij>0) then
1839 0 : LIBPAW_ALLOCATE(pawrhoij_out(irhoij)%grhoij,(ngrhoij,cplex*qphase*lmn2_size,nspden))
1840 0 : do isp=1,nspden
1841 0 : do ii=1,cplex*qphase*lmn2_size
1842 0 : pawrhoij_out(irhoij)%grhoij(1:ngrhoij,ii,isp)=buf_dp(indx_dp:indx_dp+ngrhoij-1)
1843 0 : indx_dp=indx_dp+ngrhoij
1844 : end do
1845 : end do
1846 : end if
1847 0 : if (use_rhoijres>0) then
1848 0 : LIBPAW_ALLOCATE(pawrhoij_out(irhoij)%rhoijres,(cplex*qphase*lmn2_size,nspden))
1849 0 : do isp=1,nspden
1850 : pawrhoij_out(irhoij)%rhoijres(1:cplex*qphase*lmn2_size,isp)= &
1851 0 : & buf_dp(indx_dp:indx_dp+cplex*qphase*lmn2_size-1)
1852 0 : indx_dp=indx_dp+cplex*qphase*lmn2_size
1853 : end do
1854 : end if
1855 0 : if (use_rhoij_>0) then
1856 0 : LIBPAW_ALLOCATE(pawrhoij_out(irhoij)%rhoij_,(cplex*qphase*lmn2_size,rhoij_size2))
1857 0 : do isp=1,rhoij_size2
1858 : pawrhoij_out(irhoij)%rhoij_(1:cplex*qphase*lmn2_size,isp)= &
1859 0 : & buf_dp(indx_dp:indx_dp+cplex*qphase*lmn2_size-1)
1860 0 : indx_dp=indx_dp+cplex*qphase*lmn2_size
1861 : end do
1862 : end if
1863 : end do
1864 : !Check
1865 0 : if ((indx_int/=1+buf_int_size).or.(indx_dp/=1+buf_dp_size)) then
1866 0 : msg='(3) Wrong buffer sizes !'
1867 0 : LIBPAW_BUG(msg)
1868 : end if
1869 :
1870 : !Free memory
1871 0 : LIBPAW_DEALLOCATE(buf_int)
1872 0 : LIBPAW_DEALLOCATE(buf_dp)
1873 0 : if (paral_atom) then
1874 0 : LIBPAW_POINTER_DEALLOCATE(buf_int_all)
1875 0 : LIBPAW_POINTER_DEALLOCATE(buf_dp_all)
1876 0 : LIBPAW_DEALLOCATE(count_int)
1877 0 : LIBPAW_DEALLOCATE(count_dp)
1878 0 : LIBPAW_DEALLOCATE(disp_int)
1879 0 : LIBPAW_DEALLOCATE(disp_dp)
1880 0 : LIBPAW_DEALLOCATE(atmtab)
1881 0 : call free_my_atmtab(my_atmtab,my_atmtab_allocated)
1882 : end if
1883 :
1884 0 : end subroutine pawrhoij_bcast
1885 : !!***
1886 :
1887 : !----------------------------------------------------------------------
1888 :
1889 : !!****f* m_pawrhoij/pawrhoij_redistribute
1890 : !! NAME
1891 : !! pawrhoij_redistribute
1892 : !!
1893 : !! FUNCTION
1894 : !! Redistribute an array of pawrhoij datastructures
1895 : !! Input pawrhoij is given on a MPI communicator
1896 : !! Output pawrhoij is redistributed on another MPI communicator
1897 : !!
1898 : !! INPUTS
1899 : !! mpi_comm_in= input MPI (atom) communicator
1900 : !! mpi_comm_out= output MPI (atom) communicator
1901 : !! mpi_atmtab_in= --optional-- indexes of the input pawrhoij treated by current proc
1902 : !! if not present, will be calculated in the present routine
1903 : !! mpi_atmtab_out= --optional-- indexes of the output pawrhoij treated by current proc
1904 : !! if not present, will be calculated in the present routine
1905 : !! natom= --optional-- total number of atoms
1906 : !! ----- Optional arguments used only for asynchronous communications -----
1907 : !! RecvAtomProc(:)= rank of processor from which I expect atom (in mpi_comm_in)
1908 : !! RecvAtomList(:)= indexes of atoms to be received by me
1909 : !! RecvAtomList(irecv) are the atoms I expect from RecvAtomProc(irecv)
1910 : !! SendAtomProc(:)= ranks of process destination of atom (in mpi_comm_in)
1911 : !! SendAtomList(:)= indexes of atoms to be sent by me
1912 : !! SendAtomList(isend) are the atoms sent to SendAtomProc(isend)
1913 : !!
1914 : !! OUTPUT
1915 : !! [pawrhoij_out(:)]<type(pawrhoij_type)>= --optional--
1916 : !! if present, the redistributed datastructure does not replace
1917 : !! the input one but is delivered in pawrhoij_out
1918 : !! if not present, input and output datastructure are the same.
1919 : !!
1920 : !! SIDE EFFECTS
1921 : !! pawrhoij(:)<type(pawrhoij_type)>= input (and eventually output) pawrhoij datastructures
1922 : !!
1923 : !! SOURCE
1924 :
1925 40 : subroutine pawrhoij_redistribute(pawrhoij,mpi_comm_in,mpi_comm_out,&
1926 40 : & natom,mpi_atmtab_in,mpi_atmtab_out,pawrhoij_out,&
1927 40 : & SendAtomProc,SendAtomList,RecvAtomProc,RecvAtomList)
1928 :
1929 : !Arguments ------------------------------------
1930 : !scalars
1931 : integer,intent(in) :: mpi_comm_in,mpi_comm_out
1932 : integer,optional,intent(in) :: natom
1933 : !arrays
1934 : integer,intent(in),optional,target :: mpi_atmtab_in(:),mpi_atmtab_out(:)
1935 : integer,intent(in),optional :: SendAtomProc(:),SendAtomList(:),RecvAtomProc(:),RecvAtomList(:)
1936 : type(pawrhoij_type),allocatable,intent(inout) :: pawrhoij(:)
1937 : type(pawrhoij_type),pointer,intent(out),optional :: pawrhoij_out(:)
1938 :
1939 : !Local variables-------------------------------
1940 : !scalars
1941 : integer :: algo_option,i1,iat_in,iat_out,iatom,ierr,ireq,iircv,iisend,imsg,imsg_current,imsg1
1942 : integer :: iproc_rcv,iproc_send,me_exch,mpi_comm_exch,my_natom_in,my_natom_out,my_tag,natom_tot,nb_dp,nb_int
1943 : integer :: nb_msg,nbmsg_incoming,nbrecv,nbrecvmsg,nbsend,nbsendreq,nbsent,next,npawrhoij_sent
1944 : integer :: nproc_in,nproc_out
1945 : logical :: flag,in_place,message_yet_prepared,my_atmtab_in_allocated,my_atmtab_out_allocated,paral_atom
1946 : !arrays
1947 : integer :: buf_size(3),request1(3)
1948 40 : integer,pointer :: my_atmtab_in(:),my_atmtab_out(:)
1949 40 : integer,allocatable :: atmtab_send(:),atm_indx_in(:),atm_indx_out(:),buf_int1(:),From(:),request(:)
1950 40 : integer,allocatable,target :: buf_int(:)
1951 40 : integer,pointer:: buf_ints(:)
1952 40 : logical,allocatable :: msg_pick(:)
1953 40 : real(dp),allocatable :: buf_dp1(:)
1954 40 : real(dp),allocatable,target :: buf_dp(:)
1955 40 : real(dp),pointer :: buf_dps(:)
1956 40 : type(coeffi1_type),target,allocatable :: tab_buf_int(:),tab_buf_atom(:)
1957 40 : type(coeff1_type),target,allocatable :: tab_buf_dp(:)
1958 40 : type(pawrhoij_type),pointer :: pawrhoij_out1(:)
1959 40 : type(pawrhoij_type),allocatable :: pawrhoij_all(:)
1960 :
1961 : ! *************************************************************************
1962 :
1963 : !@pawrhoij_type
1964 :
1965 40 : in_place=(.not.present(pawrhoij_out))
1966 40 : my_natom_in=size(pawrhoij)
1967 :
1968 : !If not "in_place", destroy the output datastructure
1969 40 : if (.not.in_place) then
1970 0 : if (associated(pawrhoij_out)) then
1971 0 : call pawrhoij_free(pawrhoij_out)
1972 0 : LIBPAW_DATATYPE_DEALLOCATE(pawrhoij_out)
1973 : end if
1974 : end if
1975 :
1976 : !Special sequential case
1977 40 : if (mpi_comm_in==xmpi_comm_self.and.mpi_comm_out==xmpi_comm_self) then
1978 0 : if ((.not.in_place).and.(my_natom_in>0)) then
1979 0 : LIBPAW_DATATYPE_ALLOCATE(pawrhoij_out,(my_natom_in))
1980 0 : call pawrhoij_nullify(pawrhoij_out)
1981 : call pawrhoij_copy(pawrhoij,pawrhoij_out,&
1982 0 : & keep_cplex=.false.,keep_qphase=.false.,keep_itypat=.false.,keep_nspden=.false.)
1983 : end if
1984 0 : return
1985 : end if
1986 :
1987 : !Get total natom
1988 40 : if (present(natom)) then
1989 40 : natom_tot=natom
1990 : else
1991 0 : natom_tot=my_natom_in
1992 0 : call xmpi_sum(natom_tot,mpi_comm_in,ierr)
1993 : end if
1994 :
1995 : !Select input distribution
1996 40 : if (present(mpi_atmtab_in)) then
1997 40 : my_atmtab_in => mpi_atmtab_in
1998 40 : my_atmtab_in_allocated=.false.
1999 : else
2000 : call get_my_atmtab(mpi_comm_in,my_atmtab_in,my_atmtab_in_allocated,&
2001 0 : & paral_atom,natom_tot,my_natom_in)
2002 : end if
2003 :
2004 : !Select output distribution
2005 40 : if (present(mpi_atmtab_out)) then
2006 40 : my_natom_out=size(mpi_atmtab_out)
2007 40 : my_atmtab_out => mpi_atmtab_out
2008 40 : my_atmtab_out_allocated=.false.
2009 : else
2010 : call get_my_atmtab(mpi_comm_out,my_atmtab_out,my_atmtab_out_allocated,&
2011 0 : & paral_atom,natom_tot)
2012 : end if
2013 :
2014 : !Select algo according to optional input arguments
2015 40 : algo_option=1
2016 : if (present(SendAtomProc).and.present(SendAtomList).and.&
2017 40 : & present(RecvAtomProc).and.present(RecvAtomList)) algo_option=2
2018 :
2019 :
2020 : !Brute force algorithm (allgather + scatter)
2021 : !---------------------------------------------------------
2022 : if (algo_option==1) then
2023 :
2024 0 : LIBPAW_DATATYPE_ALLOCATE(pawrhoij_all,(natom_tot))
2025 0 : call pawrhoij_nullify(pawrhoij_all)
2026 : call pawrhoij_copy(pawrhoij,pawrhoij_all,comm_atom=mpi_comm_in,mpi_atmtab=my_atmtab_in,&
2027 0 : & keep_cplex=.false.,keep_qphase=.false.,keep_itypat=.false.,keep_nspden=.false.)
2028 0 : if (in_place) then
2029 0 : call pawrhoij_free(pawrhoij)
2030 0 : LIBPAW_DATATYPE_DEALLOCATE(pawrhoij)
2031 0 : LIBPAW_DATATYPE_ALLOCATE(pawrhoij,(my_natom_out))
2032 0 : call pawrhoij_nullify(pawrhoij)
2033 : call pawrhoij_copy(pawrhoij_all,pawrhoij,comm_atom=mpi_comm_out,mpi_atmtab=my_atmtab_out,&
2034 0 : & keep_cplex=.false.,keep_qphase=.false.,keep_itypat=.false.,keep_nspden=.false.)
2035 : else
2036 0 : LIBPAW_DATATYPE_ALLOCATE(pawrhoij_out,(my_natom_out))
2037 0 : call pawrhoij_nullify(pawrhoij_out)
2038 : call pawrhoij_copy(pawrhoij_all,pawrhoij_out,comm_atom=mpi_comm_out,mpi_atmtab=my_atmtab_out,&
2039 0 : & keep_cplex=.false.,keep_qphase=.false.,keep_itypat=.false.,keep_nspden=.false.)
2040 : end if
2041 0 : call pawrhoij_free(pawrhoij_all)
2042 0 : LIBPAW_DATATYPE_DEALLOCATE(pawrhoij_all)
2043 :
2044 :
2045 : !Asynchronous algorithm (asynchronous communications)
2046 : !---------------------------------------------------------
2047 : else if (algo_option==2) then
2048 :
2049 40 : nbsend=size(SendAtomProc) ; nbrecv=size(RecvAtomProc)
2050 :
2051 40 : if (in_place) then
2052 40 : if ( my_natom_out > 0 ) then
2053 112 : LIBPAW_DATATYPE_ALLOCATE(pawrhoij_out1,(my_natom_out))
2054 28 : call pawrhoij_nullify(pawrhoij_out1)
2055 : else
2056 12 : LIBPAW_DATATYPE_ALLOCATE(pawrhoij_out1,(0))
2057 : end if
2058 : else
2059 0 : LIBPAW_DATATYPE_ALLOCATE(pawrhoij_out,(my_natom_out))
2060 0 : call pawrhoij_nullify(pawrhoij_out)
2061 0 : pawrhoij_out1=>pawrhoij_out
2062 : end if
2063 :
2064 40 : nproc_in=xmpi_comm_size(mpi_comm_in)
2065 40 : nproc_out=xmpi_comm_size(mpi_comm_out)
2066 40 : if (nproc_in<=nproc_out) mpi_comm_exch=mpi_comm_out
2067 40 : if (nproc_in>nproc_out) mpi_comm_exch=mpi_comm_in
2068 40 : me_exch=xmpi_comm_rank(mpi_comm_exch)
2069 :
2070 : ! Dimension put to the maximum to send
2071 120 : LIBPAW_ALLOCATE(atmtab_send,(nbsend))
2072 120 : LIBPAW_ALLOCATE(atm_indx_in,(natom_tot))
2073 120 : atm_indx_in=-1
2074 68 : do iatom=1,my_natom_in
2075 68 : atm_indx_in(my_atmtab_in(iatom))=iatom
2076 : end do
2077 80 : LIBPAW_ALLOCATE(atm_indx_out,(natom_tot))
2078 120 : atm_indx_out=-1
2079 68 : do iatom=1,my_natom_out
2080 68 : atm_indx_out(my_atmtab_out(iatom))=iatom
2081 : end do
2082 :
2083 148 : LIBPAW_DATATYPE_ALLOCATE(tab_buf_int,(nbsend))
2084 148 : LIBPAW_DATATYPE_ALLOCATE(tab_buf_dp,(nbsend))
2085 108 : LIBPAW_DATATYPE_ALLOCATE(tab_buf_atom,(nbsend))
2086 120 : LIBPAW_ALLOCATE(request,(3*nbsend))
2087 :
2088 : ! A send buffer in an asynchrone communication couldn't be deallocate before it has been receive
2089 40 : nbsent=0 ; ireq=0 ; iisend=0 ; nbsendreq=0 ; nb_msg=0
2090 68 : do iisend=1,nbsend
2091 28 : iproc_rcv=SendAtomProc(iisend)
2092 28 : next=-1
2093 28 : if (iisend < nbsend) next=SendAtomProc(iisend+1)
2094 68 : if (iproc_rcv /= me_exch) then
2095 8 : nbsent=nbsent+1
2096 8 : atmtab_send(nbsent)=SendAtomList(iisend) ! we groups the atoms sends to the same process
2097 8 : if (iproc_rcv /= next) then
2098 8 : if (nbsent > 0) then
2099 : ! Check if message has been yet prepared
2100 : message_yet_prepared=.false.
2101 8 : do imsg=1,nb_msg
2102 8 : if (size(tab_buf_atom(imsg)%value) /= nbsent) then
2103 : cycle
2104 : else
2105 0 : do imsg1=1,nbsent
2106 0 : if (tab_buf_atom(imsg)%value(imsg1)/=atmtab_send(imsg1)) exit
2107 0 : message_yet_prepared=.true.
2108 0 : imsg_current=imsg
2109 : end do
2110 : end if
2111 : end do
2112 : ! Create the message
2113 8 : if (.not.message_yet_prepared) then
2114 8 : nb_msg=nb_msg+1
2115 : call pawrhoij_isendreceive_fillbuffer( &
2116 8 : & pawrhoij,atmtab_send,atm_indx_in,nbsent,buf_int,nb_int,buf_dp,nb_dp)
2117 24 : LIBPAW_ALLOCATE(tab_buf_int(nb_msg)%value,(nb_int))
2118 24 : LIBPAW_ALLOCATE(tab_buf_dp(nb_msg)%value,(nb_dp))
2119 232 : tab_buf_int(nb_msg)%value(1:nb_int)=buf_int(1:nb_int)
2120 104 : tab_buf_dp(nb_msg)%value(1:nb_dp)=buf_dp(1:nb_dp)
2121 8 : LIBPAW_DEALLOCATE(buf_int)
2122 8 : LIBPAW_DEALLOCATE(buf_dp)
2123 24 : LIBPAW_ALLOCATE(tab_buf_atom(nb_msg)%value, (nbsent))
2124 16 : tab_buf_atom(nb_msg)%value(1:nbsent)=atmtab_send(1:nbsent)
2125 : imsg_current=nb_msg
2126 : end if
2127 : ! Communicate
2128 8 : buf_size(1)=size(tab_buf_int(imsg_current)%value)
2129 8 : buf_size(2)=size(tab_buf_dp(imsg_current)%value)
2130 8 : buf_size(3)=nbsent
2131 8 : buf_ints=>tab_buf_int(imsg_current)%value
2132 8 : buf_dps=>tab_buf_dp(imsg_current)%value
2133 8 : my_tag=100
2134 8 : ireq=ireq+1
2135 8 : call xmpi_isend(buf_size,iproc_rcv,my_tag,mpi_comm_exch,request(ireq),ierr)
2136 8 : my_tag=101
2137 8 : ireq=ireq+1
2138 8 : call xmpi_isend(buf_ints,iproc_rcv,my_tag,mpi_comm_exch,request(ireq),ierr)
2139 8 : my_tag=102
2140 8 : ireq=ireq+1
2141 8 : call xmpi_isend(buf_dps,iproc_rcv,my_tag,mpi_comm_exch,request(ireq),ierr)
2142 8 : nbsendreq=ireq
2143 8 : nbsent=0
2144 : end if
2145 : end if
2146 : else ! Just a renumbering, not a sending
2147 20 : iat_in=atm_indx_in(SendAtomList(iisend))
2148 20 : iat_out=atm_indx_out(my_atmtab_in(iat_in))
2149 : call pawrhoij_copy(pawrhoij(iat_in:iat_in),pawrhoij_out1(iat_out:iat_out), &
2150 20 : & keep_cplex=.false.,keep_qphase=.false.,keep_itypat=.false.,keep_nspden=.false.)
2151 20 : nbsent=0
2152 : end if
2153 : end do
2154 :
2155 120 : LIBPAW_ALLOCATE(From,(nbrecv))
2156 68 : From(:)=-1 ; nbrecvmsg=0
2157 68 : do iircv=1,nbrecv
2158 28 : iproc_send=RecvAtomProc(iircv) !receive from (RcvAtomProc is sorted by growing process)
2159 28 : next=-1
2160 28 : if (iircv < nbrecv) next=RecvAtomProc(iircv+1)
2161 68 : if (iproc_send /= me_exch .and. iproc_send/=next) then
2162 8 : nbrecvmsg=nbrecvmsg+1
2163 8 : From(nbrecvmsg)=iproc_send
2164 : end if
2165 : end do
2166 :
2167 120 : LIBPAW_ALLOCATE(msg_pick,(nbrecvmsg))
2168 48 : msg_pick=.false.
2169 : nbmsg_incoming=nbrecvmsg
2170 125 : do while (nbmsg_incoming > 0)
2171 210 : do i1=1,nbrecvmsg
2172 170 : if (.not.msg_pick(i1)) then
2173 85 : iproc_send=From(i1)
2174 : flag=.false.
2175 85 : my_tag=100
2176 85 : call xmpi_iprobe(iproc_send,my_tag,mpi_comm_exch,flag,ierr)
2177 85 : if (flag) then
2178 8 : msg_pick(i1)=.true.
2179 8 : call xmpi_irecv(buf_size,iproc_send,my_tag,mpi_comm_exch,request1(1),ierr)
2180 8 : call xmpi_wait(request1(1),ierr)
2181 8 : nb_int=buf_size(1)
2182 8 : nb_dp=buf_size(2)
2183 8 : npawrhoij_sent=buf_size(3)
2184 24 : LIBPAW_ALLOCATE(buf_int1,(nb_int))
2185 24 : LIBPAW_ALLOCATE(buf_dp1,(nb_dp))
2186 8 : my_tag=101
2187 8 : call xmpi_irecv(buf_int1,iproc_send,my_tag,mpi_comm_exch,request1(2),ierr)
2188 8 : my_tag=102
2189 8 : call xmpi_irecv(buf_dp1,iproc_send,my_tag,mpi_comm_exch,request1(3),ierr)
2190 8 : call xmpi_waitall(request1(2:3),ierr)
2191 8 : call pawrhoij_isendreceive_getbuffer(pawrhoij_out1,npawrhoij_sent,atm_indx_out,buf_int1,buf_dp1)
2192 8 : nbmsg_incoming=nbmsg_incoming-1
2193 8 : LIBPAW_DEALLOCATE(buf_int1)
2194 24 : LIBPAW_DEALLOCATE(buf_dp1)
2195 : end if
2196 : end if
2197 : end do
2198 : end do
2199 40 : LIBPAW_DEALLOCATE(msg_pick)
2200 :
2201 40 : if (in_place) then
2202 40 : call pawrhoij_free(pawrhoij)
2203 68 : LIBPAW_DATATYPE_DEALLOCATE(pawrhoij)
2204 148 : LIBPAW_DATATYPE_ALLOCATE(pawrhoij,(my_natom_out))
2205 40 : call pawrhoij_nullify(pawrhoij)
2206 : call pawrhoij_copy(pawrhoij_out1,pawrhoij, &
2207 40 : & keep_cplex=.false.,keep_qphase=.false.,keep_itypat=.false.,keep_nspden=.false.)
2208 40 : call pawrhoij_free(pawrhoij_out1)
2209 40 : LIBPAW_DATATYPE_DEALLOCATE(pawrhoij_out1)
2210 : end if
2211 :
2212 : ! Wait for deallocating arrays that all sending operations has been realized
2213 40 : if (nbsendreq > 0) then
2214 8 : call xmpi_waitall(request(1:nbsendreq),ierr)
2215 : end if
2216 :
2217 : ! Deallocate buffers
2218 48 : do i1=1,nb_msg
2219 8 : LIBPAW_DEALLOCATE(tab_buf_int(i1)%value)
2220 8 : LIBPAW_DEALLOCATE(tab_buf_dp(i1)%value)
2221 48 : LIBPAW_DEALLOCATE(tab_buf_atom(i1)%value)
2222 : end do
2223 68 : LIBPAW_DATATYPE_DEALLOCATE(tab_buf_int)
2224 68 : LIBPAW_DATATYPE_DEALLOCATE(tab_buf_dp)
2225 68 : LIBPAW_DATATYPE_DEALLOCATE(tab_buf_atom)
2226 40 : LIBPAW_DEALLOCATE(From)
2227 40 : LIBPAW_DEALLOCATE(request)
2228 40 : LIBPAW_DEALLOCATE(atmtab_send)
2229 40 : LIBPAW_DEALLOCATE(atm_indx_in)
2230 40 : LIBPAW_DEALLOCATE(atm_indx_out)
2231 :
2232 : end if !algo_option
2233 :
2234 : !Eventually release temporary pointers
2235 40 : call free_my_atmtab(my_atmtab_in,my_atmtab_in_allocated)
2236 40 : call free_my_atmtab(my_atmtab_out,my_atmtab_out_allocated)
2237 :
2238 80 : end subroutine pawrhoij_redistribute
2239 : !!***
2240 :
2241 : !----------------------------------------------------------------------
2242 :
2243 : !!****f* m_pawrhoij/pawrhoij_io
2244 : !! NAME
2245 : !! pawrhoij_io
2246 : !!
2247 : !! FUNCTION
2248 : !! IO method for pawrhoij datastructures.
2249 : !!
2250 : !! INPUTS
2251 : !! unitfi=Unit number for IO file or netcdf file handler (already opened in the caller).
2252 : !! nsppol_in=Number of independent spin polarizations. Only used for reading.
2253 : !! nspinor_in=Number of spinorial components. Only used for reading.
2254 : !! nspden_in=Number of spin-density components. only used for reading.
2255 : !! nlmn_type(ntypat)= Number of (l,m,n) elements for the paw basis for each type of atom. Only used for reading.
2256 : !! typat(natom) =Type of each atom.
2257 : !! headform=Format of the abinit header (only used for reading as we need to know how to read
2258 : !! the data. Writing is always done using the latest headform.
2259 : !! rdwr_mode(len=*)=String defining the IO mode. Possible values (not case sensitive):
2260 : !! "W"= For writing to unitfi
2261 : !! "R"= For reading from unitfi
2262 : !! "E"= For echoing.
2263 : !! "D"= for debug
2264 : !! [form(len=*)]= String defining the file format. Defaults to Fortran binary mode i.e., "unformatted"
2265 : !! Other possible values are (case insensitive):
2266 : !! "formatted"=For IO on a file open in formatted mode.
2267 : !! "netcdf"=For IO on a netcdf file.
2268 : !! [natinc]=Defines the increment in the loop over natom used for echoing the pawrhoij(natom) datastructures.
2269 : !! If not specified, only the first and the last atom will be printed.
2270 : !!
2271 : !! SIDE EFFECTS
2272 : !! pawrhoij(:)<type(pawrhoij_type)>= rhoij datastructure.
2273 : !! if rdwr_mode="W", it will be written on unit unitfi using the file format given by form.
2274 : !! if rdwr_mode="R", pawrhoij will be read and initialized from unit unitfi that has been
2275 : !! opened with form=form.
2276 : !! if rdwr_mode="E", the routines only echoes the content of the structure.
2277 : !!
2278 : !! SOURCE
2279 :
2280 5421 : subroutine pawrhoij_io(pawrhoij,unitfi,nsppol_in,nspinor_in,nspden_in,nlmn_type,typat,&
2281 : & headform,rdwr_mode,form,natinc,mpi_atmtab)
2282 :
2283 : !Arguments ------------------------------------
2284 : !scalars
2285 : integer,intent(in) :: unitfi,headform,nspden_in,nspinor_in,nsppol_in
2286 : integer,optional,intent(in) :: natinc
2287 : character(len=*),intent(in) :: rdwr_mode
2288 : character(len=*),optional,intent(in) :: form
2289 : integer, intent(in), optional, pointer :: mpi_atmtab(:)
2290 : !arrays
2291 : integer,intent(in) :: typat(:),nlmn_type(:)
2292 : type(pawrhoij_type),intent(inout),target :: pawrhoij(:)
2293 :
2294 : !Local variables-------------------------------
2295 : !scalars
2296 : integer,parameter :: fort_formatted=1,fort_binary=2,netcdf_io=3
2297 : integer :: cplex,qphase,i1,i2,iatom,iatom_tot,natom,ispden,bsize,ii,jj,lmn2_size
2298 : integer :: nselect,my_cplex,my_cplex_eff,my_qphase,my_natinc,my_natom,my_nspden,ngrhoijmx,size_rhoij2
2299 : integer :: iomode,ncid,natom_id,cplex_id,qphase_id,nspden_id,nsel56_id
2300 : integer :: buffer_id,ibuffer_id,ncerr,bsize_id,bufsize_id
2301 : integer :: iq0,itypat,lmn_size
2302 : logical :: paral_atom
2303 : character(len=500) :: msg
2304 : !arrays
2305 5421 : integer,allocatable :: ibuffer(:),nsel44(:,:),nsel56(:)
2306 5421 : integer,pointer :: my_atmtab(:)
2307 5421 : real(dp), allocatable :: buffer(:)
2308 5421 : real(dp),pointer :: rhoij_tmp(:)
2309 :
2310 : ! *************************************************************************
2311 :
2312 5421 : my_natom=SIZE(pawrhoij);if (my_natom==0) return
2313 5421 : my_nspden=nspden_in
2314 5421 : natom=size(typat)
2315 5421 : paral_atom=(my_natom/=natom)
2316 5421 : if (present(mpi_atmtab)) then
2317 0 : if (.not.associated(mpi_atmtab)) then
2318 0 : msg='mpi_atmtab not associated (pawrhoij_io)'
2319 0 : LIBPAW_BUG(msg)
2320 : end if
2321 0 : my_atmtab=>mpi_atmtab
2322 5421 : else if (my_natom/=natom) then
2323 0 : msg='my_natom /=natom, mpi_atmtab should be in argument (pawrhoij_io)'
2324 0 : LIBPAW_BUG(msg)
2325 : end if
2326 :
2327 5421 : iomode = fort_binary
2328 5421 : if (PRESENT(form)) then
2329 10526 : select case (libpaw_to_upper(form))
2330 : case ("FORMATTED")
2331 5263 : iomode = fort_formatted
2332 : case ("NETCDF")
2333 5263 : iomode = netcdf_io
2334 : case default
2335 5263 : LIBPAW_ERROR("Wrong form: "//trim(form))
2336 : end select
2337 : end if
2338 :
2339 : #ifndef LIBPAW_HAVE_NETCDF
2340 : if (iomode == netcdf_io) then
2341 : LIBPAW_ERROR("iomode == netcdf_io but netcdf library is missing.")
2342 : end if
2343 : #endif
2344 5421 : ncid = unitfi
2345 :
2346 2964 : select case (rdwr_mode(1:1))
2347 :
2348 : case ("R","r") ! Reading the Rhoij tab
2349 :
2350 2964 : if ((headform>=44).and.(headform<56)) then
2351 0 : LIBPAW_ALLOCATE(nsel44,(nspden_in,natom))
2352 0 : if (iomode == fort_binary) then
2353 0 : read(unitfi ) ((nsel44(ispden,iatom),ispden=1,nspden_in),iatom=1,natom)
2354 0 : else if (iomode == fort_formatted) then
2355 0 : read(unitfi,*) ((nsel44(ispden,iatom),ispden=1,nspden_in),iatom=1,natom)
2356 : #ifdef LIBPAW_HAVE_NETCDF
2357 : else if (iomode == netcdf_io) then
2358 0 : LIBPAW_ERROR("header in 44-56 not compatible with Netcdf")
2359 : #endif
2360 : end if
2361 0 : call pawrhoij_alloc(pawrhoij,1,nspden_in,nspinor_in,nsppol_in,typat,lmnsize=nlmn_type)
2362 0 : do iatom=1,natom
2363 0 : pawrhoij(iatom)%nrhoijsel=nsel44(1,iatom)
2364 : end do
2365 0 : bsize=sum(nsel44)
2366 0 : LIBPAW_ALLOCATE(ibuffer,(bsize))
2367 0 : LIBPAW_ALLOCATE(buffer,(bsize))
2368 0 : if (iomode == fort_binary) then
2369 0 : read(unitfi ) ibuffer(:),buffer(:)
2370 0 : else if (iomode == fort_formatted) then
2371 0 : read(unitfi,*) ibuffer(:),buffer(:)
2372 : end if
2373 : ii=0
2374 0 : do iatom=1,natom
2375 0 : nselect=nsel44(1,iatom)
2376 0 : pawrhoij(iatom)%rhoijselect(:)=0
2377 0 : pawrhoij(iatom)%rhoijselect(1:nselect)=ibuffer(ii+1:ii+nselect)
2378 0 : do ispden=1,nspden_in
2379 0 : pawrhoij(iatom)%rhoijp(1:nselect,ispden)=buffer(ii+1:ii+nselect)
2380 0 : ii=ii+nselect
2381 : end do
2382 : end do
2383 0 : LIBPAW_DEALLOCATE(ibuffer)
2384 0 : LIBPAW_DEALLOCATE(buffer)
2385 0 : LIBPAW_DEALLOCATE(nsel44)
2386 2964 : else if (headform>=56) then
2387 8892 : LIBPAW_ALLOCATE(nsel56,(natom))
2388 2964 : my_cplex=1;my_nspden=1;my_qphase=1
2389 2964 : if (headform==56) then
2390 0 : if (iomode == fort_binary) then
2391 0 : read(unitfi ) (nsel56(iatom),iatom=1,natom),my_cplex
2392 0 : else if (iomode == fort_formatted) then
2393 0 : read(unitfi,*) (nsel56(iatom),iatom=1,natom),my_cplex
2394 : #ifdef LIBPAW_HAVE_NETCDF
2395 : else if (iomode == netcdf_io) then
2396 0 : NCF_CHECK(nf90_inq_dimid(ncid, "pawrhoij_cplex", cplex_id))
2397 0 : NCF_CHECK(nf90_inquire_dimension(ncid, cplex_id, len=my_cplex))
2398 :
2399 0 : NCF_CHECK(nf90_inq_varid(ncid, "rhoijsel_atoms", nsel56_id))
2400 0 : NCF_CHECK(nf90_get_var(ncid, nsel56_id, nsel56))
2401 : #endif
2402 : end if
2403 : else
2404 2964 : if (iomode == fort_binary) then
2405 6 : read(unitfi,err=10,end=10) (nsel56(iatom),iatom=1,natom),my_cplex,my_nspden,my_qphase
2406 : 10 continue
2407 2958 : else if (iomode == fort_formatted) then
2408 0 : read(unitfi,fmt=*,err=11,end=11) (nsel56(iatom),iatom=1,natom),my_cplex,my_nspden,my_qphase
2409 : 11 continue
2410 : #ifdef LIBPAW_HAVE_NETCDF
2411 : else if (iomode == netcdf_io) then
2412 2958 : NCF_CHECK(nf90_inq_dimid(ncid, "pawrhoij_cplex", cplex_id))
2413 2958 : NCF_CHECK(nf90_inquire_dimension(ncid, cplex_id, len=my_cplex))
2414 2958 : NCF_CHECK(nf90_inq_dimid(ncid, "pawrhoij_nspden", nspden_id))
2415 2958 : NCF_CHECK(nf90_inquire_dimension(ncid, nspden_id, len=my_nspden))
2416 2958 : NCF_CHECK(nf90_inq_varid(ncid, "nrhoijsel_atoms", nsel56_id))
2417 2958 : NCF_CHECK(nf90_get_var(ncid, nsel56_id, nsel56))
2418 2958 : if (nf90_inq_dimid(ncid, "pawrhoij_qphase", qphase_id)==NF90_NOERR) then
2419 2958 : NCF_CHECK(nf90_inquire_dimension(ncid, qphase_id, len=my_qphase))
2420 : else
2421 : my_qphase=1
2422 : end if
2423 : #endif
2424 : end if
2425 : end if
2426 : call pawrhoij_alloc(pawrhoij,my_cplex,my_nspden,nspinor_in,nsppol_in,typat,&
2427 2964 : & lmnsize=nlmn_type,qphase=my_qphase)
2428 10100 : do iatom=1,natom
2429 10100 : pawrhoij(iatom)%nrhoijsel=nsel56(iatom)
2430 : end do
2431 10100 : bsize=sum(nsel56)
2432 8892 : LIBPAW_ALLOCATE(ibuffer,(bsize))
2433 8892 : LIBPAW_ALLOCATE(buffer,(bsize*my_nspden*my_cplex*my_qphase))
2434 2964 : if (iomode == fort_binary) then
2435 6 : read(unitfi ) ibuffer(:),buffer(:)
2436 2958 : else if (iomode == fort_formatted) then
2437 0 : read(unitfi,*) ibuffer(:),buffer(:)
2438 : #ifdef LIBPAW_HAVE_NETCDF
2439 : else if (iomode == netcdf_io) then
2440 2958 : if (bsize > 0) then
2441 2092 : NCF_CHECK(nf90_inq_varid(ncid, "rhoijselect_atoms", ibuffer_id))
2442 2092 : NCF_CHECK(nf90_get_var(ncid, ibuffer_id, ibuffer))
2443 2092 : NCF_CHECK(nf90_inq_varid(ncid, "rhoijp_atoms", buffer_id))
2444 2092 : NCF_CHECK(nf90_get_var(ncid, buffer_id, buffer))
2445 : end if
2446 : #endif
2447 : end if
2448 : ii=0;jj=0
2449 10100 : do iatom=1,natom
2450 7136 : nselect=nsel56(iatom)
2451 316609 : pawrhoij(iatom)%rhoijselect(:)=0
2452 147154 : pawrhoij(iatom)%rhoijselect(1:nselect)=ibuffer(ii+1:ii+nselect)
2453 7136 : ii=ii+nselect
2454 20885 : do ispden=1,my_nspden
2455 : pawrhoij(iatom)%rhoijp(1:my_cplex*nselect,ispden)= &
2456 281279 : buffer(jj+1:jj+my_cplex*nselect)
2457 10785 : jj=jj+my_cplex*nselect
2458 17921 : if (my_qphase==2) then
2459 0 : itypat=typat(iatom)
2460 0 : lmn_size=nlmn_type(itypat)
2461 0 : lmn2_size=lmn_size*(lmn_size+1)/2
2462 0 : iq0 = my_cplex*lmn2_size
2463 : pawrhoij(iatom)%rhoijp(iq0+1:iq0+my_cplex*nselect,ispden)= &
2464 0 : buffer(jj+1:jj+my_cplex*nselect)
2465 : jj=jj+my_cplex*nselect
2466 : end if
2467 : end do
2468 : end do
2469 2964 : LIBPAW_DEALLOCATE(ibuffer)
2470 2964 : LIBPAW_DEALLOCATE(buffer)
2471 2964 : LIBPAW_DEALLOCATE(nsel56)
2472 : end if
2473 :
2474 : case ("W","w") ! Writing the Rhoij tab (latest format is used)
2475 :
2476 7365 : LIBPAW_ALLOCATE(nsel56,(natom))
2477 2455 : my_cplex =pawrhoij(1)%cplex_rhoij
2478 2455 : my_nspden=pawrhoij(1)%nspden
2479 2455 : my_qphase=pawrhoij(1)%qphase
2480 8520 : do iatom=1,natom
2481 8520 : nsel56(iatom)=pawrhoij(iatom)%nrhoijsel
2482 : end do
2483 8520 : bsize=sum(nsel56)
2484 :
2485 2455 : if (iomode == fort_binary) then
2486 150 : write(unitfi ) (nsel56(iatom),iatom=1,natom),my_cplex,my_nspden,my_qphase
2487 2305 : else if (iomode == fort_formatted) then
2488 0 : write(unitfi,*) (nsel56(iatom),iatom=1,natom),my_cplex,my_nspden,my_qphase
2489 : #ifdef LIBPAW_HAVE_NETCDF
2490 : else if (iomode == netcdf_io) then
2491 2305 : ncerr = nf90_redef(ncid)
2492 :
2493 : ! Define dimensions.
2494 2305 : ncerr = nf90_inq_dimid(ncid, "number_of_atoms", natom_id)
2495 2305 : if (ncerr /= nf90_noerr) then
2496 0 : NCF_CHECK(nf90_def_dim(ncid, "number_of_atoms", natom, natom_id))
2497 : end if
2498 2305 : ncerr = nf90_inq_varid(ncid, "nrhoijsel_atoms", nsel56_id)
2499 2305 : if (ncerr /= nf90_noerr) then
2500 2304 : NCF_CHECK(nf90_def_var(ncid, "nrhoijsel_atoms", NF90_INT, natom_id, nsel56_id))
2501 : end if
2502 2305 : ncerr = nf90_inq_dimid(ncid, "pawrhoij_cplex", cplex_id)
2503 2305 : if (ncerr /= nf90_noerr) then
2504 2304 : NCF_CHECK(nf90_def_dim(ncid, "pawrhoij_cplex", my_cplex, cplex_id))
2505 : end if
2506 2305 : ncerr = nf90_inq_dimid(ncid, "pawrhoij_nspden", nspden_id)
2507 2305 : if (ncerr /= nf90_noerr) then
2508 2304 : NCF_CHECK(nf90_def_dim(ncid, "pawrhoij_nspden", my_nspden, nspden_id))
2509 : end if
2510 2305 : ncerr = nf90_inq_dimid(ncid, "pawrhoij_qphase", qphase_id)
2511 2305 : if (ncerr /= nf90_noerr) then
2512 2304 : NCF_CHECK(nf90_def_dim(ncid, "pawrhoij_qphase", my_qphase, qphase_id))
2513 : end if
2514 2305 : if (bsize > 0) then
2515 1999 : ncerr = nf90_inq_dimid(ncid, "rhoijselect_atoms_dim", bsize_id)
2516 1999 : if (ncerr /= nf90_noerr) then
2517 1998 : NCF_CHECK(nf90_def_dim(ncid, "rhoijselect_atoms_dim", bsize, bsize_id))
2518 : end if
2519 1999 : ncerr = nf90_inq_dimid(ncid, "rhoijp_atoms_dim", bufsize_id)
2520 1999 : if (ncerr /= nf90_noerr) then
2521 1998 : NCF_CHECK(nf90_def_dim(ncid, "rhoijp_atoms_dim", bsize*my_nspden*my_qphase*my_cplex, bufsize_id))
2522 : end if
2523 1999 : ncerr = nf90_inq_varid(ncid, "rhoijselect_atoms", ibuffer_id)
2524 1999 : if (ncerr /= nf90_noerr) then
2525 1998 : NCF_CHECK(nf90_def_var(ncid, "rhoijselect_atoms", NF90_INT, bsize_id, ibuffer_id))
2526 : end if
2527 1999 : ncerr = nf90_inq_varid(ncid, "rhoijp_atoms", buffer_id)
2528 1999 : if (ncerr /= nf90_noerr) then
2529 1998 : NCF_CHECK(nf90_def_var(ncid, "rhoijp_atoms", NF90_DOUBLE, bufsize_id, buffer_id))
2530 : end if
2531 : else
2532 : ! This happens in v5[40] and bsize == 0 corresponds to NC_UNLIMITED
2533 306 : LIBPAW_COMMENT("All rhoij entries are zero. No netcdf entry produced")
2534 : end if
2535 :
2536 : ! Write nsel56
2537 2305 : NCF_CHECK(nf90_enddef(ncid))
2538 2305 : NCF_CHECK(nf90_put_var(ncid, nsel56_id, nsel56))
2539 : #endif
2540 : end if
2541 :
2542 7365 : LIBPAW_ALLOCATE(ibuffer,(bsize))
2543 7365 : LIBPAW_ALLOCATE(buffer,(bsize*my_nspden*my_cplex*my_qphase))
2544 2455 : ii=0;jj=0
2545 8520 : do iatom=1,natom
2546 6065 : nselect=nsel56(iatom)
2547 228092 : ibuffer(ii+1:ii+nselect)=pawrhoij(iatom)%rhoijselect(1:nselect)
2548 6065 : ii=ii+nselect
2549 17205 : do ispden=1,my_nspden
2550 : buffer(jj+1:jj+my_cplex*nselect)=&
2551 516795 : pawrhoij(iatom)%rhoijp(1:my_cplex*nselect,ispden)
2552 8685 : jj=jj+my_cplex*nselect
2553 14750 : if (my_qphase==2) then
2554 152 : iq0 = my_cplex*pawrhoij(iatom)%lmn2_size
2555 : buffer(jj+1:jj+my_cplex*nselect)=&
2556 8040 : pawrhoij(iatom)%rhoijp(iq0+1:iq0+my_cplex*nselect,ispden)
2557 : jj=jj+my_cplex*nselect
2558 : end if
2559 : end do
2560 : end do
2561 2455 : if (iomode == fort_binary) then
2562 150 : write(unitfi ) ibuffer(:),buffer(:)
2563 2305 : else if (iomode == fort_formatted) then
2564 0 : write(unitfi,*) ibuffer(:),buffer(:)
2565 : #ifdef LIBPAW_HAVE_NETCDF
2566 : else if (iomode == netcdf_io) then
2567 2305 : if (bsize > 0) then
2568 1999 : NCF_CHECK(nf90_put_var(ncid, ibuffer_id, ibuffer))
2569 1999 : NCF_CHECK(nf90_put_var(ncid, buffer_id, buffer))
2570 : end if
2571 : #endif
2572 : end if
2573 2455 : LIBPAW_DEALLOCATE(ibuffer)
2574 2455 : LIBPAW_DEALLOCATE(buffer)
2575 2455 : LIBPAW_DEALLOCATE(nsel56)
2576 :
2577 : case ("E","e") ! Echoing the Rhoij tab
2578 :
2579 2 : my_natinc=1; if(natom>1) my_natinc=natom-1
2580 2 : my_qphase=pawrhoij(1)%qphase
2581 : nselect=maxval(pawrhoij(:)%nrhoijsel)
2582 2 : if (PRESENT(natinc)) my_natinc = natinc ! user-defined increment.
2583 2 : LIBPAW_ALLOCATE(ibuffer,(0))
2584 6 : nselect=maxval(pawrhoij(:)%nrhoijsel)
2585 2 : if (my_qphase==2) then
2586 0 : LIBPAW_POINTER_ALLOCATE(rhoij_tmp,(2*nselect))
2587 : end if
2588 6 : do iatom=1,my_natom,my_natinc
2589 4 : iatom_tot=iatom;if(paral_atom)iatom_tot=my_atmtab(iatom)
2590 4 : my_cplex =pawrhoij(iatom)%cplex_rhoij
2591 4 : my_nspden=pawrhoij(iatom)%nspden
2592 4 : my_qphase=pawrhoij(iatom)%qphase
2593 4 : nselect=pawrhoij(iatom)%nrhoijsel
2594 10 : do ispden=1,pawrhoij(iatom)%nspden
2595 4 : if (my_qphase==1) then
2596 4 : my_cplex_eff=my_cplex
2597 4 : rhoij_tmp => pawrhoij(iatom)%rhoijp(1:my_cplex*nselect,ispden)
2598 : else
2599 0 : my_cplex_eff=2
2600 0 : rhoij_tmp=zero
2601 0 : jj=my_cplex*pawrhoij(iatom)%lmn2_size
2602 0 : if (my_cplex==1) then
2603 0 : do ii=1,nselect
2604 0 : rhoij_tmp(2*ii-1)=pawrhoij(iatom)%rhoijp(ii,ispden)
2605 0 : rhoij_tmp(2*ii )=pawrhoij(iatom)%rhoijp(jj+ii,ispden)
2606 : end do
2607 : else
2608 0 : do ii=1,nselect
2609 : rhoij_tmp(2*ii-1)=pawrhoij(iatom)%rhoijp(2*ii-1,ispden) &
2610 0 : & -pawrhoij(iatom)%rhoijp(jj+2*ii ,ispden)
2611 : rhoij_tmp(2*ii )=pawrhoij(iatom)%rhoijp(2*ii ,ispden) &
2612 0 : & +pawrhoij(iatom)%rhoijp(jj+2*ii-1,ispden)
2613 : end do
2614 : end if
2615 : end if
2616 4 : write(unitfi, '(a,i4,a,i1,a)' ) ' rhoij(',iatom_tot,',',ispden,')= (max 12 non-zero components will be written)'
2617 : call pawio_print_ij(unitfi,rhoij_tmp,nselect,my_cplex_eff,&
2618 : & pawrhoij(iatom)%lmn_size,-1,ibuffer,1,0,&
2619 : & pawrhoij(iatom)%rhoijselect,-1.d0,1,&
2620 8 : & opt_sym=2,mode_paral='PERS')
2621 : end do ! end nspden do
2622 : end do ! end iatom do
2623 2 : LIBPAW_DEALLOCATE(ibuffer)
2624 2 : if (my_qphase==2) then
2625 0 : LIBPAW_POINTER_DEALLOCATE(rhoij_tmp)
2626 : end if
2627 :
2628 : case ("D","d") ! Debug
2629 :
2630 0 : write(unitfi,'(a,i4)' ) 'size pawmkrhoij , natom = ' , my_natom
2631 0 : my_natinc=1; if(natom>1) my_natinc=natom-1
2632 0 : if (PRESENT(natinc)) my_natinc = natinc ! user-defined increment.
2633 0 : LIBPAW_ALLOCATE(ibuffer,(0))
2634 0 : do iatom=1,my_natom,my_natinc
2635 0 : iatom_tot=iatom;if(paral_atom) iatom_tot=my_atmtab(iatom)
2636 0 : if (iatom_tot/=1) cycle
2637 0 : write(unitfi,'(a,i4,a)' ) ' ******* rhoij (Atom # ',iatom_tot,' ********)'
2638 0 : write(unitfi,'(a,i4,a,i4)' ) 'cplex_rhoij=',pawrhoij(iatom)%cplex_rhoij, ' nselect=', pawrhoij(iatom)%nrhoijsel
2639 0 : write(unitfi,'(a,i4,a,i4)' ) 'nspden=',pawrhoij(iatom)%nspden, ' lmn2size=',pawrhoij(iatom)%lmn2_size
2640 0 : write(unitfi,'(a,i4,a,i4)' ) 'lmnmix=',pawrhoij(iatom)%lmnmix_sz, ' ngrhoij=',pawrhoij(iatom)%ngrhoij
2641 0 : write(unitfi,'(a,i4,a,i4)' ) 'use_rhoijres=',pawrhoij(iatom)%use_rhoijres, &
2642 0 : & 'use_rhoij_=',pawrhoij(iatom)%use_rhoij_
2643 0 : write(unitfi,'(a,i4,a,i4)' ) 'itypat=',pawrhoij(iatom)%itypat, ' lmn_size=',pawrhoij(iatom)%lmn_size
2644 0 : write(unitfi,'(a,i4,a,i4)' ) 'nsppol=',pawrhoij(iatom)%nsppol, ' nspinor=',pawrhoij(iatom)%nspinor
2645 0 : write(unitfi,'(a,i4)' ) 'qphase=',pawrhoij(iatom)%qphase
2646 0 : cplex=pawrhoij(iatom)%cplex_rhoij
2647 0 : qphase=pawrhoij(iatom)%qphase
2648 0 : lmn2_size=pawrhoij(iatom)%lmn2_size
2649 0 : do i2=1,pawrhoij(iatom)%nrhoijsel
2650 0 : write(unitfi,'(a,i4,a,i4,a,i4,a,f9.5)') 'rhoijselect(,',i2,')=',&
2651 0 : & pawrhoij(iatom)%rhoijselect(i2)
2652 : end do
2653 0 : if (pawrhoij(iatom)%ngrhoij>0) then
2654 0 : ngrhoijmx=2; if (pawrhoij(iatom)%ngrhoij<ngrhoijmx) ngrhoijmx=pawrhoij(iatom)%ngrhoij
2655 0 : do ispden=1,pawrhoij(iatom)%nspden
2656 0 : do i1=ngrhoijmx,ngrhoijmx
2657 0 : do ii=1,qphase
2658 0 : do i2=(ii-1)*cplex*lmn2_size+cplex*lmn2_size,(ii-1)*cplex*lmn2_size+cplex*lmn2_size
2659 0 : write(unitfi,'(a,i4,a,i4,a,i4,a,f9.5)') ' grhoij(,',i1,',',i2,',',ispden,')=',&
2660 0 : & pawrhoij(iatom)%grhoij(i1,i2,ispden)
2661 : end do
2662 : end do
2663 : end do
2664 : end do
2665 0 : call libpaw_flush(unitfi)
2666 : end if
2667 0 : if (pawrhoij(iatom)%use_rhoijres>0) then
2668 0 : do ispden=1,pawrhoij(iatom)%nspden
2669 0 : do ii=1,qphase
2670 0 : do i2=(ii-1)*cplex*lmn2_size+cplex*lmn2_size,(ii-1)*cplex*lmn2_size+cplex*lmn2_size
2671 0 : write(unitfi,'(a,i4,a,i4,a,f9.5)') ' rhoijres(,',i2,',ispden=',ispden,')=',&
2672 0 : & pawrhoij(iatom)%rhoijres(i2,ispden)
2673 : end do
2674 : end do
2675 : end do
2676 0 : call libpaw_flush(unitfi)
2677 : end if
2678 0 : if (pawrhoij(iatom)%nrhoijsel>0) then
2679 0 : do ispden=1,pawrhoij(iatom)%nspden
2680 0 : do ii=1,qphase
2681 0 : do i2=(ii-1)*cplex*lmn2_size+cplex*pawrhoij(iatom)%nrhoijsel, &
2682 0 : & (ii-1)*cplex*lmn2_size+cplex*pawrhoij(iatom)%nrhoijsel
2683 0 : write(unitfi,'(a,i4,a,i4,a,f9.5)') ' rhoijp!(nrhoijselec=,',i2,',ispden=',ispden,')=',&
2684 0 : & pawrhoij(iatom)%rhoijp(i2,ispden)
2685 : end do
2686 : end do
2687 : end do
2688 0 : call libpaw_flush(unitfi)
2689 : end if
2690 0 : if (pawrhoij(iatom)%use_rhoij_>0) then
2691 0 : size_rhoij2=size(pawrhoij(iatom)%rhoij_,2)
2692 0 : do ispden=1,size_rhoij2
2693 0 : do ii=1,qphase
2694 0 : do i2=(ii-1)*cplex*lmn2_size+cplex*lmn2_size,(ii-1)*cplex*lmn2_size+cplex*lmn2_size
2695 0 : write(unitfi,'(a,i4,a,i4,a,f9.5)') ' rhoij_(,',i2,',ispden=',ispden,')=',&
2696 0 : & pawrhoij(iatom)%rhoij_(i2,ispden)
2697 : end do
2698 : end do
2699 : end do
2700 : end if
2701 0 : call libpaw_flush(unitfi)
2702 0 : if (pawrhoij(iatom)%lmnmix_sz>0) then
2703 0 : write(unitfi,'(a)') 'kpawmix '
2704 0 : write(unitfi,'(i4,i4,i4,i4)') pawrhoij(iatom)%kpawmix(:)
2705 : end if
2706 0 : call libpaw_flush(unitfi)
2707 : end do
2708 :
2709 : case default
2710 0 : msg='Wrong rdwr_mode'//TRIM(rdwr_mode)
2711 5421 : LIBPAW_ERROR(msg)
2712 :
2713 : end select
2714 :
2715 10842 : end subroutine pawrhoij_io
2716 : !!***
2717 :
2718 : !----------------------------------------------------------------------
2719 :
2720 : !!****f* m_pawrhoij/pawrhoij_unpack
2721 : !! NAME
2722 : !! pawrhoij_unpack
2723 : !!
2724 : !! FUNCTION
2725 : !! Unpack the values store in rhoijp copying them to the rhoij_ array.
2726 : !!
2727 : !! SIDE EFFECTS
2728 : !! rhoij(:) <pawrhoij_type)>= input/output datastructure
2729 : !! * In output the rhoij_ array is filled with the values stored in the packed array rhoijp.
2730 : !! * If use_rhoij_/=1, rhoij_ is allocated and the corresponding flag is set to 1.
2731 : !!
2732 : !! SOURCE
2733 :
2734 0 : subroutine pawrhoij_unpack(rhoij)
2735 :
2736 : !Arguments ------------------------------------
2737 : !scalars
2738 : !arrays
2739 : type(pawrhoij_type),intent(inout) :: rhoij(:)
2740 :
2741 : !Local variables-------------------------------
2742 : integer :: cplex_rhoij,natom,lmn2_size,nsp2,qphase
2743 : integer :: i0,iat,iphase,isel,isppol,klmn
2744 :
2745 : ! *************************************************************************
2746 :
2747 0 : natom = SIZE(rhoij) ; if (natom==0) return
2748 :
2749 0 : do iat=1,natom
2750 :
2751 0 : lmn2_size =rhoij(iat)%lmn2_size
2752 0 : cplex_rhoij = rhoij(iat)%cplex_rhoij
2753 0 : qphase = rhoij(iat)%qphase
2754 0 : nsp2 = rhoij(iat)%nsppol;if (rhoij(iat)%nspden==4) nsp2=4
2755 :
2756 0 : if (rhoij(iat)%use_rhoij_/=1) then ! Have to allocate rhoij
2757 0 : LIBPAW_ALLOCATE(rhoij(iat)%rhoij_,(cplex_rhoij*qphase*lmn2_size,nsp2))
2758 0 : rhoij(iat)%use_rhoij_=1
2759 : end if
2760 0 : rhoij(iat)%rhoij_ = zero
2761 :
2762 0 : do isppol=1,nsp2
2763 0 : do iphase=1,qphase
2764 0 : i0=(iphase-1)*lmn2_size*cplex_rhoij
2765 0 : if (cplex_rhoij==1) then
2766 0 : do isel=1,rhoij(iat)%nrhoijsel ! Looping over non-zero ij elements.
2767 0 : klmn = rhoij(iat)%rhoijselect(isel)
2768 0 : rhoij(iat)%rhoij_(i0+klmn,isppol) = rhoij(iat)%rhoijp(i0+isel,isppol)
2769 : end do
2770 : else
2771 0 : do isel=1,rhoij(iat)%nrhoijsel ! Looping over non-zero ij elements.
2772 0 : klmn = rhoij(iat)%rhoijselect(isel)
2773 : rhoij(iat)%rhoij_(i0+2*klmn-1:i0+2*klmn,isppol) = &
2774 0 : & rhoij(iat)%rhoijp(i0+2*isel-1:i0+2*isel,isppol)
2775 : end do
2776 : end if
2777 : end do
2778 : end do
2779 :
2780 : end do ! natom
2781 :
2782 : end subroutine pawrhoij_unpack
2783 : !!***
2784 :
2785 : !----------------------------------------------------------------------
2786 :
2787 : !!****f* m_pawrhoij/pawrhoij_init_unpacked
2788 : !! NAME
2789 : !! pawrhoij_init_unpacked
2790 : !!
2791 : !! FUNCTION
2792 : !! Initialize field of rhoij datastructure for unpacked values (pawrhoij%rhoij_ array)
2793 : !!
2794 : !! SIDE EFFECTS
2795 : !! rhoij(:) <pawrhoij_type)>= input/output datastructure
2796 : !! * In output the rhoij_ array is allocated
2797 : !!
2798 : !! SOURCE
2799 :
2800 19434 : subroutine pawrhoij_init_unpacked(rhoij)
2801 :
2802 : !Arguments ------------------------------------
2803 : !scalars
2804 : !arrays
2805 : type(pawrhoij_type),intent(inout) :: rhoij(:)
2806 :
2807 : !Local variables-------------------------------
2808 : integer :: cplex_rhoij,iat,lmn2_size,nrhoij,nsp2,qphase
2809 :
2810 : ! *************************************************************************
2811 :
2812 19434 : nrhoij=SIZE(rhoij);if (nrhoij==0) return
2813 :
2814 72875 : do iat=1,nrhoij
2815 :
2816 53441 : lmn2_size =rhoij(iat)%lmn2_size
2817 53441 : cplex_rhoij = rhoij(iat)%cplex_rhoij
2818 53441 : qphase = rhoij(iat)%qphase
2819 53441 : nsp2 = rhoij(iat)%nsppol;if (rhoij(iat)%nspden==4) nsp2=4
2820 :
2821 53441 : if (allocated(rhoij(iat)%rhoij_)) then
2822 0 : LIBPAW_DEALLOCATE(rhoij(iat)%rhoij_)
2823 : end if
2824 213764 : LIBPAW_ALLOCATE(rhoij(iat)%rhoij_,(cplex_rhoij*qphase*lmn2_size,nsp2))
2825 53441 : rhoij(iat)%use_rhoij_=1
2826 6041839 : rhoij(iat)%rhoij_=zero
2827 :
2828 : end do
2829 :
2830 : end subroutine pawrhoij_init_unpacked
2831 : !!***
2832 :
2833 : !----------------------------------------------------------------------
2834 :
2835 : !!****f* m_pawrhoij/pawrhoij_free_unpacked
2836 : !! NAME
2837 : !! pawrhoij_free_unpacked
2838 : !!
2839 : !! FUNCTION
2840 : !! Destroy field of rhoij datastructure for unpacked values (pawrhoij%rhoij_ array)
2841 : !!
2842 : !! SIDE EFFECTS
2843 : !! rhoij(:) <pawrhoij_type)>= input/output datastructure
2844 : !! * In output the rhoij_ array is deallocated
2845 : !!
2846 : !! SOURCE
2847 :
2848 22510 : subroutine pawrhoij_free_unpacked(rhoij)
2849 :
2850 : !Arguments ------------------------------------
2851 : !scalars
2852 : !arrays
2853 : type(pawrhoij_type),intent(inout) :: rhoij(:)
2854 :
2855 : !Local variables-------------------------------
2856 : integer :: iat,nrhoij
2857 :
2858 : ! *************************************************************************
2859 :
2860 22510 : nrhoij=SIZE(rhoij);if (nrhoij==0) return
2861 :
2862 89919 : do iat=1,nrhoij
2863 :
2864 67409 : if (allocated(rhoij(iat)%rhoij_)) then
2865 66949 : LIBPAW_DEALLOCATE(rhoij(iat)%rhoij_)
2866 : end if
2867 89919 : rhoij(iat)%use_rhoij_=0
2868 :
2869 : end do
2870 :
2871 : end subroutine pawrhoij_free_unpacked
2872 : !!***
2873 :
2874 : !----------------------------------------------------------------------
2875 :
2876 : !!****f* m_pawrhoij/pawrhoij_mpisum_unpacked_1D
2877 : !! NAME
2878 : !! pawrhoij_mpisum_unpacked_1D
2879 : !!
2880 : !! FUNCTION
2881 : !! Build the MPI sum of the unsymmetrized PAW rhoij_ (augmentation occupancies)
2882 : !! Remember:for each atom, rho_ij=Sum_{n,k} {occ(n,k)*<Cnk|p_i><p_j|Cnk>}
2883 : !! Target: 1D array of pawrhoij datastructures
2884 : !!
2885 : !! INPUTS
2886 : !! comm1=MPI communicator. Data will be MPI summed inside comm1
2887 : !! [comm2]=second MPI communicator. If present, rhoij_ will be
2888 : !! MPI summed inside comm2 after the collective sum in comm1.
2889 : !!
2890 : !! SIDE EFFECTS
2891 : !! pawrhoij(:) <type(pawrhoij_type)>= paw rhoij occupancies and related data
2892 : !! Input: the data calculateed by this processor.
2893 : !! Output: the final MPI sum over comm1 and comm2.
2894 : !!
2895 : !! SOURCE
2896 :
2897 19566 : subroutine pawrhoij_mpisum_unpacked_1D(pawrhoij,comm1,comm2)
2898 :
2899 : !Arguments ---------------------------------------------
2900 : !scalars
2901 : integer,intent(in) :: comm1
2902 : integer,optional,intent(in) :: comm2
2903 : !arrays
2904 : type(pawrhoij_type),intent(inout) :: pawrhoij(:)
2905 :
2906 : !Local variables ---------------------------------------
2907 : !scalars
2908 : integer :: bufdim,iatom,ierr,isppol,jdim,nsp2,natom
2909 : integer :: nproc1,nproc2
2910 : !character(len=500) :: msg
2911 : !arrays
2912 19566 : integer,allocatable :: dimlmn(:)
2913 19566 : real(dp),allocatable :: buffer(:)
2914 :
2915 : !************************************************************************
2916 :
2917 19566 : natom=SIZE(pawrhoij);if (natom==0) return
2918 :
2919 19566 : nproc1 = xmpi_comm_size(comm1)
2920 19566 : nproc2=1; if (PRESENT(comm2)) nproc2 = xmpi_comm_size(comm2)
2921 19566 : if (nproc1==1.and.nproc2==1) RETURN
2922 :
2923 : !Fill the MPI buffer from the local rhoij_
2924 22692 : LIBPAW_ALLOCATE(dimlmn,(natom))
2925 : dimlmn(1:natom)=pawrhoij(1:natom)%cplex_rhoij &
2926 : & *pawrhoij(1:natom)%qphase &
2927 27580 : & *pawrhoij(1:natom)%lmn2_size
2928 7564 : nsp2=pawrhoij(1)%nsppol; if (pawrhoij(1)%nspden==4) nsp2=4
2929 27580 : bufdim=sum(dimlmn)*nsp2
2930 22692 : LIBPAW_ALLOCATE(buffer,(bufdim))
2931 7564 : jdim=0
2932 27580 : do iatom=1,natom
2933 53144 : do isppol=1,nsp2
2934 1877108 : buffer(jdim+1:jdim+dimlmn(iatom))=pawrhoij(iatom)%rhoij_(:,isppol)
2935 45580 : jdim=jdim+dimlmn(iatom)
2936 : end do
2937 : end do
2938 :
2939 : !Build sum of pawrhoij%rhoij_
2940 7564 : call xmpi_sum(buffer,comm1,ierr) ! Sum over the first communicator.
2941 7564 : if (PRESENT(comm2)) then
2942 5132 : call xmpi_sum(buffer,comm2,ierr) ! Sum over the second communicator.
2943 : end if
2944 :
2945 : !Unpack the MPI packet filling rhoij_
2946 : jdim=0
2947 27580 : do iatom=1,natom
2948 53144 : do isppol=1,nsp2
2949 1877108 : pawrhoij(iatom)%rhoij_(:,isppol)=buffer(jdim+1:jdim+dimlmn(iatom))
2950 45580 : jdim=jdim+dimlmn(iatom)
2951 : end do
2952 : end do
2953 :
2954 7564 : LIBPAW_DEALLOCATE(buffer)
2955 7564 : LIBPAW_DEALLOCATE(dimlmn)
2956 :
2957 7564 : end subroutine pawrhoij_mpisum_unpacked_1D
2958 : !!***
2959 :
2960 : !----------------------------------------------------------------------
2961 :
2962 : !!****f* m_pawrhoij/pawrhoij_mpisum_unpacked_2D
2963 : !! NAME
2964 : !! pawrhoij_mpisum_unpacked_2D
2965 : !!
2966 : !! FUNCTION
2967 : !! Build the MPI sum of the unsymmetrized PAW rhoij_ (augmentation occupancies)
2968 : !! Remember:for each atom, rho_ij=Sum_{n,k} {occ(n,k)*<Cnk|p_i><p_j|Cnk>}
2969 : !! Target: 2D array of pawrhoij datastructures
2970 : !!
2971 : !! INPUTS
2972 : !! comm1=MPI communicator. Data will be MPI summed inside comm1
2973 : !! [comm2]=second MPI communicator. If present, rhoij_ will be
2974 : !! MPI summed inside comm2 after the collective sum in comm1.
2975 : !!
2976 : !! SIDE EFFECTS
2977 : !! pawrhoij(:,:) <type(pawrhoij_type)>= paw rhoij occupancies and related data
2978 : !! Input: the data calculateed by this processor.
2979 : !! Output: the final MPI sum over comm1 and comm2.
2980 : !!
2981 : !! SOURCE
2982 :
2983 942 : subroutine pawrhoij_mpisum_unpacked_2D(pawrhoij,comm1,comm2)
2984 :
2985 : !Arguments ---------------------------------------------
2986 : !scalars
2987 : integer,intent(in) :: comm1
2988 : integer,optional,intent(in) :: comm2
2989 : !arrays
2990 : type(pawrhoij_type),intent(inout) :: pawrhoij(:,:)
2991 :
2992 : !Local variables ---------------------------------------
2993 : !scalars
2994 : integer :: bufdim,iatom,ierr,irhoij,isppol,jdim,nsp2,natom,nrhoij
2995 : integer :: nproc1,nproc2
2996 : !character(len=500) :: msg
2997 : !arrays
2998 942 : integer,allocatable :: dimlmn(:,:)
2999 942 : real(dp),allocatable :: buffer(:)
3000 :
3001 : !************************************************************************
3002 :
3003 942 : natom =SIZE(pawrhoij,1);if (natom ==0) return
3004 942 : nrhoij=SIZE(pawrhoij,2);if (nrhoij==0) return
3005 :
3006 942 : nproc1 = xmpi_comm_size(comm1)
3007 942 : nproc2=1; if (PRESENT(comm2)) nproc2 = xmpi_comm_size(comm2)
3008 942 : if (nproc1==1.and.nproc2==1) RETURN
3009 :
3010 : !Fill the MPI buffer from the local rhoij_
3011 2192 : LIBPAW_ALLOCATE(dimlmn,(natom,nrhoij))
3012 : dimlmn(1:natom,1:nrhoij)=pawrhoij(1:natom,1:nrhoij)%cplex_rhoij &
3013 : & *pawrhoij(1:natom,1:nrhoij)%qphase &
3014 10052 : & *pawrhoij(1:natom,1:nrhoij)%lmn2_size
3015 548 : nsp2=pawrhoij(1,1)%nsppol; if (pawrhoij(1,1)%nspden==4) nsp2=4
3016 10052 : bufdim=sum(dimlmn)*nsp2
3017 1644 : LIBPAW_ALLOCATE(buffer,(bufdim))
3018 548 : jdim=0
3019 2192 : do irhoij=1,nrhoij
3020 10052 : do iatom=1,natom
3021 17364 : do isppol=1,nsp2
3022 297732 : buffer(jdim+1:jdim+dimlmn(iatom,irhoij))=pawrhoij(iatom,irhoij)%rhoij_(:,isppol)
3023 15720 : jdim=jdim+dimlmn(iatom,irhoij)
3024 : end do
3025 : end do
3026 : end do
3027 :
3028 : !Build sum of pawrhoij%rhoij_
3029 548 : call xmpi_sum(buffer,comm1,ierr) ! Sum over the first communicator.
3030 548 : if (PRESENT(comm2)) then
3031 0 : call xmpi_sum(buffer,comm2,ierr) ! Sum over the second communicator.
3032 : end if
3033 :
3034 : !Unpack the MPI packet filling rhoij_
3035 : jdim=0
3036 2192 : do irhoij=1,nrhoij
3037 10052 : do iatom=1,natom
3038 17364 : do isppol=1,nsp2
3039 297732 : pawrhoij(iatom,irhoij)%rhoij_(:,isppol)=buffer(jdim+1:jdim+dimlmn(iatom,irhoij))
3040 15720 : jdim=jdim+dimlmn(iatom,irhoij)
3041 : end do
3042 : end do
3043 : end do
3044 :
3045 548 : LIBPAW_DEALLOCATE(buffer)
3046 548 : LIBPAW_DEALLOCATE(dimlmn)
3047 :
3048 548 : end subroutine pawrhoij_mpisum_unpacked_2D
3049 : !!***
3050 :
3051 : !----------------------------------------------------------------------
3052 :
3053 : !!****f* m_pawrhoij/pawrhoij_filter
3054 : !! NAME
3055 : !! pawrhoij_filter
3056 : !!
3057 : !! FUNCTION
3058 : !! Filter a "rhoij" array (PAW on-site occupancies),
3059 : !! i.e. select only the non-zero elements.
3060 : !!
3061 : !! INPUT
3062 : !! cplex=1 if Rhoij is real, 2 if Rhoij is complex
3063 : !! qphase=2 if rhoij has a exp(iqR) phase, 1 if not
3064 : !! lmn2_size=dimension of rhoij=number of (i,j) pairs with i<=j
3065 : !! nspden=number of rhoij spin components
3066 : !! [rhoij_input(cplex*qphase*lmn2_size,nspden)]= -- optional argument--
3067 : !! input values for rhoij (including zero values)
3068 : !! If this argument is not provided, that the input values from rhoij()
3069 : !!
3070 : !! OUTPUT
3071 : !! nselect=number of non-zero values of rhoij
3072 : !! rhoijselect(lmn2_size)=contains the indices of the selected (i,j) pairs
3073 : !! rhoijselect(nselect+1:lmn2_size)=0
3074 : !!
3075 : !! SIDE EFFECTS
3076 : !! rhoij(cplex*qphase*lmn2_size,nspden)=
3077 : !! * Input: the array is filled with all rhoij values (only if rhoij_input is not present)
3078 : !! * Output: the nselect first elements contain the non-zero rhoij values,
3079 : !! next value are irrelevant
3080 : !!
3081 : !! SOURCE
3082 :
3083 171988 : subroutine pawrhoij_filter(rhoij,rhoijselect,nselect,cplex,qphase,lmn2_size,nspden, &
3084 55690 : & rhoij_input) ! optional argument
3085 :
3086 : !Arguments ------------------------------------
3087 : !scalars
3088 : integer,intent(in) :: lmn2_size,cplex,qphase,nspden
3089 : integer,intent(out) :: nselect
3090 : !arrays
3091 : integer,intent(out) :: rhoijselect(lmn2_size)
3092 : real(dp),intent(inout),target :: rhoij(cplex*qphase*lmn2_size,nspden)
3093 : real(dp),intent(in),optional,target :: rhoij_input(cplex*qphase*lmn2_size,nspden)
3094 :
3095 : !Local variables-------------------------------
3096 : !scalars
3097 : real(dp),parameter :: tol_rhoij=tol10
3098 : integer :: isp,klmn,klmn1,klmn2,nsel1,nsel2
3099 : !arrays
3100 85994 : real(dp),pointer :: rhoij_in(:,:)
3101 :
3102 : ! *************************************************************************
3103 :
3104 85994 : nselect=0
3105 4612686 : rhoijselect(:)=0
3106 :
3107 85994 : if (present(rhoij_input)) then
3108 55690 : rhoij_in => rhoij_input
3109 : else
3110 30304 : rhoij_in => rhoij
3111 : end if
3112 :
3113 : !Treat each cplex/qphase case separately
3114 :
3115 85994 : if (cplex==1) then
3116 :
3117 80855 : if (qphase==1) then
3118 :
3119 4062192 : do klmn=1,lmn2_size
3120 6555472 : if (any(abs(rhoij_in(klmn,:))>tol_rhoij)) then
3121 2261081 : nselect=nselect+1
3122 2261081 : rhoijselect(nselect)=klmn
3123 5509060 : do isp=1,nspden
3124 5509060 : rhoij(nselect,isp)=rhoij_in(klmn,isp)
3125 : end do
3126 : end if
3127 : end do
3128 :
3129 1380 : else if (qphase==2) then
3130 :
3131 55920 : do klmn=1,lmn2_size
3132 54540 : klmn2=klmn+lmn2_size
3133 68577 : if (any(abs(rhoij_in(klmn,:))>tol_rhoij).or. &
3134 1380 : & any(abs(rhoij_in(klmn2,:))>tol_rhoij)) then
3135 50252 : nselect=nselect+1 ; nsel2=nselect+lmn2_size
3136 50252 : rhoijselect(nselect)=klmn
3137 106307 : do isp=1,nspden
3138 56055 : rhoij(nselect,isp)=rhoij_in(klmn ,isp)
3139 106307 : rhoij(nsel2 ,isp)=rhoij_in(klmn2,isp)
3140 : end do
3141 : end if
3142 : end do
3143 :
3144 : end if
3145 :
3146 : else ! cplex=2
3147 :
3148 5139 : if (qphase==1) then
3149 494574 : do klmn=1,lmn2_size
3150 489435 : klmn1=2*klmn
3151 1211667 : if (any(abs(rhoij_in(klmn1-1:klmn1,:))>tol_rhoij)) then
3152 462622 : nselect=nselect+1 ; nsel1=2*nselect
3153 462622 : rhoijselect(nselect)=klmn
3154 2245818 : do isp=1,nspden
3155 1783196 : rhoij(nsel1-1,isp)=rhoij_in(klmn1-1,isp)
3156 2245818 : rhoij(nsel1 ,isp)=rhoij_in(klmn1 ,isp)
3157 : end do
3158 : end if
3159 : end do
3160 :
3161 0 : else if (qphase==2) then
3162 :
3163 0 : do klmn=1,lmn2_size
3164 0 : klmn1=2*klmn ; klmn2=klmn1+lmn2_size
3165 0 : if (any(abs(rhoij_in(klmn1-1:klmn1,:))>tol_rhoij).or. &
3166 0 : & any(abs(rhoij_in(klmn2-1:klmn2,:))>tol_rhoij)) then
3167 0 : nselect=nselect+1 ; nsel1=2*nselect ; nsel2=nsel1+lmn2_size
3168 0 : rhoijselect(nselect)=klmn
3169 0 : do isp=1,nspden
3170 0 : rhoij(nsel1-1,isp)=rhoij_in(klmn1-1,isp)
3171 0 : rhoij(nsel1 ,isp)=rhoij_in(klmn1 ,isp)
3172 0 : rhoij(nsel2-1,isp)=rhoij_in(klmn2-1,isp)
3173 0 : rhoij(nsel2 ,isp)=rhoij_in(klmn2 ,isp)
3174 : end do
3175 : end if
3176 : end do
3177 :
3178 : end if
3179 : endif
3180 :
3181 141684 : end subroutine pawrhoij_filter
3182 : !!***
3183 :
3184 : !----------------------------------------------------------------------
3185 :
3186 : !!****f* m_pawrhoij/pawrhoij_inquire_dim
3187 : !! NAME
3188 : !! pawrhoij_inquire_dim
3189 : !!
3190 : !! FUNCTION
3191 : !! Compute the values f the dimensions (cplex_rhoij, qphase, nspden) for a pawrhoij datastructure.
3192 : !! These ones depend on the parameters of the calculation
3193 : !!
3194 : !! INPUTS
3195 : !! [cplex]= flag controlling the use of a exp(iqR) phase. 1=no phase, 2=phase
3196 : !! [cpxocc]= 2 if rhoij is required to be imaginary
3197 : !! [nspden]= number of spin-density components
3198 : !! [qpt(3)]= q-vector, if any
3199 : !! [spnorb]= flag: 1 if spin-orbit coupling is activated
3200 : !!
3201 : !! OUTPUT
3202 : !! [cplex_rhoij] = value of cplex_rhoij associated to pawrhoij
3203 : !! [qphase_rhoij]= value of qphase associated to pawrhoij
3204 : !! [nspden_rhoij]= value of nspden associated to pawrhoij
3205 : !!
3206 : !! SOURCE
3207 :
3208 13412 : subroutine pawrhoij_inquire_dim(cplex,cpxocc,nspden,qpt,spnorb, &
3209 : & cplex_rhoij,qphase_rhoij,nspden_rhoij)
3210 :
3211 : !Arguments ---------------------------------------------
3212 : !scalars
3213 : integer,optional,intent(in) :: cplex,cpxocc,nspden,spnorb
3214 : integer,optional,intent(out) :: cplex_rhoij,qphase_rhoij,nspden_rhoij
3215 : !arrays
3216 : real(dp),optional,intent(in) :: qpt(3)
3217 :
3218 : !Local variables ---------------------------------------
3219 : character(len=100) :: msg
3220 :
3221 : !************************************************************************
3222 :
3223 : !cplex_rhoij
3224 13412 : if (present(cplex_rhoij)) then
3225 13412 : cplex_rhoij=1
3226 13412 : if (present(cpxocc)) cplex_rhoij=max(cplex_rhoij,cpxocc)
3227 : end if
3228 :
3229 : !qphase_rhoij
3230 13412 : if (present(qphase_rhoij)) then
3231 8232 : qphase_rhoij=1
3232 8232 : if (present(cplex).and.present(qpt)) then
3233 0 : msg='only one argument cplex or qpt should be passed!'
3234 0 : LIBPAW_BUG(msg)
3235 : end if
3236 8638 : if (present(cplex)) qphase_rhoij=merge(1,2,cplex==1)
3237 8232 : if (present(qpt)) then
3238 8718 : if (any(abs(qpt(:))>tol8)) qphase_rhoij=2
3239 : end if
3240 : end if
3241 :
3242 : !nspden_rhoij
3243 13412 : if (present(nspden_rhoij)) then
3244 13412 : nspden_rhoij=1
3245 13412 : if (present(nspden)) nspden_rhoij=nspden
3246 14528 : if (present(spnorb)) nspden_rhoij=merge(nspden_rhoij,4,spnorb<=0)
3247 : end if
3248 :
3249 13412 : end subroutine pawrhoij_inquire_dim
3250 : !!***
3251 :
3252 : !----------------------------------------------------------------------
3253 :
3254 : !!****f* m_pawdij/pawrhoij_print_rhoij
3255 : !! NAME
3256 : !! pawrhoij_print_rhoij
3257 : !!
3258 : !! FUNCTION
3259 : !! Print out the content of a Rho_ij matrix (occupation matrix) in a suitable format
3260 : !!
3261 : !! INPUTS
3262 : !! rhoij(cplex*lmn2_size,nspden)= input matrix to be printed
3263 : !! cplex=1 if Rhoij is real, 2 if Rhoij is complex
3264 : !! qphase=2 if rhoij has a exp(iqR) phase, 1 if not
3265 : !! iatom=current atom
3266 : !! natom=total number of atoms in the system
3267 : !! [opt_prtvol]= >=0 if up to 12 components of _ij matrix have to be printed
3268 : !! <0 if all components of ij_ matrix have to be printed (optional)
3269 : !! [mode_paral]= parallel printing mode (optional, default='COLL')
3270 : ! [rhoijselect(lmn2_size)]=Indirect array selecting the non-zero elements of rhoij
3271 : !! [test_value]=(real number) if positive, print a warning when the magnitude of Dij is greater (optional)
3272 : !! [l_only]=if >=0 only parts of rhoij corresponding to li=l_only are printed (optional);
3273 : !! Needs indlmn(:,:) optional argument.
3274 : !! [title_msg]=message to print as title (optional)
3275 : !! [unit]=the unit number for output (optional)
3276 : !! [indlmn(6,lmn_size)]= array giving l,m,n,lm,ln,s
3277 : !!
3278 : !! OUTPUT
3279 : !! (Only writing)
3280 : !!
3281 : !! NOTES
3282 : !!
3283 : !! SOURCE
3284 :
3285 21870 : subroutine pawrhoij_print_rhoij(rhoij,cplex,qphase,iatom,natom,&
3286 21870 : & rhoijselect,test_value,title_msg,unit,opt_prtvol,&
3287 21870 : & l_only,indlmn,mode_paral) ! Optional arguments
3288 :
3289 : !Arguments ------------------------------------
3290 : !scalars
3291 : integer,intent(in) :: cplex,qphase,iatom,natom
3292 : integer,optional,intent(in) :: opt_prtvol,l_only,unit
3293 : real(dp),intent(in),optional :: test_value
3294 : character(len=4),optional,intent(in) :: mode_paral
3295 : character(len=100),optional,intent(in) :: title_msg
3296 : !arrays
3297 : integer,optional,intent(in) :: indlmn(:,:)
3298 : integer,optional,intent(in),target :: rhoijselect(:)
3299 : real(dp),intent(in),target :: rhoij(:,:)
3300 :
3301 : !Local variables-------------------------------
3302 : !scalars
3303 : character(len=8),parameter :: dspin(6)=(/"up ","down ","dens (n)","magn (x)","magn (y)","magn (z)"/)
3304 : integer :: irhoij,kk,my_cplex,my_lmn_size,my_lmn2_size,my_l_only,my_nspden
3305 : integer :: my_opt_pack,my_opt_sym,my_prtvol,my_unt,nrhoijsel,rhoij_size
3306 : real(dp) :: my_test_value,test_value_eff
3307 : character(len=4) :: my_mode
3308 : character(len=2000) :: msg
3309 : !arrays
3310 : integer,target :: idum(0)
3311 21870 : integer,pointer :: l_index(:),my_rhoijselect(:)
3312 21870 : real(dp),pointer :: rhoij_(:)
3313 :
3314 : ! *************************************************************************
3315 :
3316 : !Optional arguments
3317 21870 : my_unt =std_out ; if (PRESENT(unit )) my_unt =unit
3318 21870 : my_mode ='COLL' ; if (PRESENT(mode_paral)) my_mode =mode_paral
3319 21870 : my_prtvol=1 ; if (PRESENT(opt_prtvol)) my_prtvol=opt_prtvol
3320 21870 : my_test_value=-one; if (PRESENT(test_value)) my_test_value=test_value
3321 21870 : my_l_only=-1 ; if (PRESENT(l_only)) my_l_only=l_only
3322 :
3323 21870 : if (my_l_only>=0.and.(.not.present(indlmn))) then
3324 0 : msg='pawrhoij_print_rhoij: l_only>=0 and indlmn not present!'
3325 0 : LIBPAW_BUG(msg)
3326 : end if
3327 :
3328 : !Title
3329 21870 : if (present(title_msg)) then
3330 0 : if (trim(title_msg)/='') then
3331 0 : write(msg, '(2a)') ch10,trim(title_msg)
3332 0 : call wrtout(my_unt,msg,my_mode)
3333 : end if
3334 : end if
3335 :
3336 : !Inits
3337 21870 : my_nspden=size(rhoij,2)
3338 21870 : my_lmn2_size=size(rhoij,1)/cplex/qphase
3339 21870 : my_lmn_size=int(dsqrt(two*dble(my_lmn2_size)))
3340 21870 : my_cplex=merge(cplex,2,qphase==1)
3341 21870 : my_opt_sym=2
3342 :
3343 : !Packed storage
3344 21870 : my_opt_pack=0
3345 21870 : rhoij_size=my_lmn2_size
3346 21870 : my_rhoijselect => idum
3347 21870 : if (PRESENT(rhoijselect)) then
3348 1000797 : nrhoijsel=count(rhoijselect(:)>0)
3349 13162 : if (nrhoijsel>0) then
3350 12472 : my_opt_pack=1
3351 12472 : rhoij_size=nrhoijsel
3352 12472 : my_rhoijselect => rhoijselect(1:nrhoijsel)
3353 : end if
3354 : end if
3355 :
3356 21870 : if (my_l_only<0) then
3357 21440 : l_index => idum
3358 : else
3359 1290 : LIBPAW_POINTER_ALLOCATE(l_index,(my_lmn_size))
3360 8036 : do kk=1,my_lmn_size
3361 8036 : l_index(kk)=indlmn(1,kk)
3362 : end do
3363 : end if
3364 :
3365 21870 : if (qphase==2) then
3366 1800 : LIBPAW_DATATYPE_ALLOCATE(rhoij_,(2*rhoij_size))
3367 : end if
3368 :
3369 : ! === Loop over Rho_ij components ===
3370 52294 : do irhoij=1,my_nspden
3371 :
3372 : !Rebuild rhoij according to qphase
3373 30424 : if (qphase==1) then
3374 29800 : rhoij_ => rhoij(1:cplex*rhoij_size,irhoij)
3375 : else
3376 624 : if (cplex==1) then
3377 27518 : do kk=1,rhoij_size
3378 26894 : rhoij_(2*kk-1)=rhoij(kk,irhoij)
3379 27518 : rhoij_(2*kk )=rhoij(kk+my_lmn2_size,irhoij)
3380 : end do
3381 : else
3382 0 : do kk=1,rhoij_size
3383 0 : rhoij_(2*kk-1)=rhoij(2*kk-1,irhoij)-rhoij(2*kk +2*my_lmn2_size,irhoij)
3384 0 : rhoij_(2*kk )=rhoij(2*kk ,irhoij)+rhoij(2*kk-1+2*my_lmn2_size,irhoij)
3385 : end do
3386 : end if
3387 : end if
3388 :
3389 : !Subtitle
3390 30424 : if (natom>1.or.my_nspden>1) then
3391 29426 : if (my_l_only<0) then
3392 28532 : if (my_nspden==1) write(msg,'(a,i3)') ' Atom #',iatom
3393 28532 : if (my_nspden==2) write(msg,'(a,i3,a,i1)')' Atom #',iatom,' - Spin component ',irhoij
3394 28532 : if (my_nspden==4) write(msg,'(a,i3,2a)') ' Atom #',iatom,' - Component ',trim(dspin(irhoij+2*(my_nspden/4)))
3395 : else
3396 908 : if (my_nspden==1) write(msg,'(a,i3,a,i1,a)') ' Atom #',iatom,&
3397 28 : & ' - L=',my_l_only,' ONLY'
3398 1494 : if (my_nspden==2) write(msg,'(a,i3,a,i1,a,i1)')' Atom #',iatom,&
3399 1200 : & ' - L=',my_l_only,' ONLY - Spin component ',irhoij
3400 1174 : if (my_nspden==4) write(msg,'(a,i3,a,i1,3a)') ' Atom #',iatom,&
3401 560 : & ' - L=',my_l_only,' ONLY - Component ',trim(dspin(irhoij+2*(my_nspden/4)))
3402 : end if
3403 29426 : call wrtout(my_unt,msg,my_mode)
3404 998 : else if (my_l_only>=0) then
3405 46 : write(msg,'(a,i1,a)') ' L=',my_l_only,' ONLY'
3406 46 : call wrtout(my_unt,msg,my_mode)
3407 : end if
3408 :
3409 : !Printing
3410 30424 : test_value_eff=-one;if(my_test_value>zero.and.irhoij==1) test_value_eff=my_test_value
3411 : call pawio_print_ij(my_unt,rhoij_,rhoij_size,my_cplex,my_lmn_size,my_l_only,l_index,my_opt_pack,&
3412 : & my_prtvol,my_rhoijselect,test_value_eff,1,opt_sym=my_opt_sym,&
3413 52294 : & mode_paral=my_mode,force_print=.true.)
3414 :
3415 : end do !irhoij
3416 :
3417 21870 : if (qphase==2) then
3418 600 : LIBPAW_DATATYPE_DEALLOCATE(rhoij_)
3419 : end if
3420 21870 : if (my_l_only>=0) then
3421 430 : LIBPAW_POINTER_DEALLOCATE(l_index)
3422 : end if
3423 :
3424 21870 : end subroutine pawrhoij_print_rhoij
3425 : !!***
3426 :
3427 : !----------------------------------------------------------------------
3428 :
3429 : !!****f* m_pawrhoij/pawrhoij_symrhoij
3430 : !! NAME
3431 : !! pawrhoij_symrhoij
3432 : !!
3433 : !! FUNCTION
3434 : !! Symmetrize rhoij quantities (augmentation occupancies) and/or gradients
3435 : !! Compute also rhoij residuals (new-old values of rhoij and gradients)
3436 : !!
3437 : !! INPUTS
3438 : !! choice=select then type of rhoij gradients to symmetrize.
3439 : !! choice=1 => no gradient
3440 : !! choice=2 => gradient with respect to atomic position(s)
3441 : !! =3 => a gradient with respect to strain(s)
3442 : !! =4 => 2nd gradient with respect to atomic position(s)
3443 : !! =23=> a gradient with respect to atm. pos. and strain(s)
3444 : !! =24=> 1st and 2nd gradient with respect to atomic position(s)
3445 : !! gprimd(3,3)=dimensional primitive translations for reciprocal space(bohr^-1).
3446 : !! indsym(4,nsym,natom)=indirect indexing array for atom labels
3447 : !! ipert=index of perturbation if pawrhoij is a pertubed rhoij
3448 : !! no meaning for ground-state calculations (should be 0)
3449 : !! [mpi_atmtab(:)]=--optional-- indexes of the atoms treated by current proc
3450 : !! [comm_atom]=--optional-- MPI communicator over atoms
3451 : !! natom=number of atoms in cell
3452 : !! nsym=number of symmetry elements in space group
3453 : !! ntypat=number of types of atoms in unit cell.
3454 : !! optrhoij= 1 if rhoij quantities have to be symmetrized
3455 : !! pawrhoij_unsym(:) <type(pawrhoij_type)>= datastructure containing PAW rhoij occupancies
3456 : !! (and related data) non symmetrized in an unpacked storage (pawrhoij_unsym%rhoij_)
3457 : !! Contains eventually unsymmetrized rhoij gradients (grhoij)
3458 : !! pawang <type(pawang_type)>=angular mesh discretization and related data
3459 : !! pawprtvol=control print volume and debugging output for PAW
3460 : !! Note: if pawprtvol=-10001, nothing is printed out
3461 : !! pawtab(ntypat) <type(pawtab_type)>=paw tabulated starting data
3462 : !! [qphon(3)]=--optional-- (RF calculations only) - wavevector of the phonon
3463 : !! rprimd(3,3)=real space primitive translations.
3464 : !! symafm(nsym)=(anti)ferromagnetic part of symmetry operations
3465 : !! symrec(3,3,nsym)=symmetries of group in terms of operations on
3466 : !! reciprocal space primitive translations
3467 : !! typat(natom)=type for each atom
3468 : !! [use_zeromag]=--optional-- .TRUE. if rhoij "magnetization" is enforced to be zero
3469 : !! Applies only when nspden_rhoij=4 (note: only the real part is set to zero)
3470 : !!
3471 : !! OUTPUT
3472 : !!
3473 : !! SIDE EFFECTS
3474 : !! pawrhoij(natom) <type(pawrhoij_type)>= datastructure containing PAW rhoij occupancies
3475 : !! (and related data) SYMMETRIZED in a packed storage (pawrhoij%rhoijp).
3476 : !! if (optrhoij==1)
3477 : !! pawrhoij(natom)%nrhoijsel=number of non-zero values of rhoij
3478 : !! pawrhoij(natom)%rhoijp(:,:)=symetrized paw rhoij quantities in PACKED STORAGE (only non-zero values)
3479 : !! pawrhoij(natom)%rhoijres(:,:)=paw rhoij quantities residuals (new values - old values)
3480 : !! pawrhoij(natom)%rhoijselect(:)=select the non-zero values of rhoij!!
3481 : !! if (pawrhoij(:)%ngrhoij>0) (equivalent to choice>1)
3482 : !! pawrhoij(natom)%grhoij(:,:)=symetrized gradients of rhoij
3483 : !!
3484 : !! NOTES
3485 : !! pawrhoij and pawrhoij_unsym can be identical (refer to the same pawrhoij datastructure).
3486 : !! They should be different only if pawrhoij is distributed over atomic sites
3487 : !! (in that case pawrhoij_unsym should not be distributed over atomic sites).
3488 : !!
3489 : !! SOURCE
3490 :
3491 22583 : subroutine pawrhoij_symrhoij(pawrhoij,pawrhoij_unsym,choice,gprimd,indsym,ipert,natom,nsym,&
3492 22583 : & ntypat,optrhoij,pawang,pawprtvol,pawtab,rprimd,symafm,symrec,typat, &
3493 22583 : & mpi_atmtab,comm_atom,qphon,use_zeromag) ! optional arguments (parallelism)
3494 :
3495 : !Arguments ---------------------------------------------
3496 : !scalars
3497 : integer,intent(in) :: choice,ipert,natom,nsym,ntypat,optrhoij,pawprtvol
3498 : integer,optional,intent(in) :: comm_atom
3499 : logical,optional,intent(in) :: use_zeromag
3500 : type(pawang_type),intent(in) :: pawang
3501 : !arrays
3502 : integer,intent(in) :: indsym(4,nsym,natom)
3503 : integer,optional,target,intent(in) :: mpi_atmtab(:)
3504 : integer,intent(in) :: symafm(nsym),symrec(3,3,nsym),typat(natom)
3505 : real(dp),intent(in) :: gprimd(3,3),rprimd(3,3)
3506 : real(dp),intent(in),optional :: qphon(3)
3507 : type(pawrhoij_type),intent(inout) :: pawrhoij(:)
3508 : type(pawrhoij_type),target,intent(inout) :: pawrhoij_unsym(:)
3509 : type(pawtab_type),target,intent(in) :: pawtab(ntypat)
3510 :
3511 : !Local variables ---------------------------------------
3512 : !scalars
3513 : integer :: at_indx,cplex_eff,cplex_rhoij,iafm,iatm,iatom,il,il0,ilmn,iln,iln0,ilpm,indexi
3514 : integer :: indexii,indexj,indexjj,indexjj0,indexk,indexkc,indexkc_q,iplex,iq,iq0
3515 : integer :: irhoij,irot,ishift2,ishift3,ishift4,ispden,itypat
3516 : integer :: j0lmn,jl,jl0,jlmn,jln,jln0,jlpm,jrhoij,jspden,klmn,klmn1,klmn1q,kspden
3517 : integer :: lmn_size,lmn2_size,mi,mj,my_comm_atom,mu,mua,mub,mushift
3518 : integer :: natinc,ngrhoij,nrhoij,nrhoij1,nrhoij_unsym
3519 : integer :: nselect,nu,nushift,qphase,sz1,sz2
3520 : real(dp) :: det
3521 : logical,parameter :: afm_noncoll=.true. ! TRUE if antiferro symmetries are used with non-collinear magnetism
3522 : logical :: use_zeromag_
3523 : real(dp) :: arg,factafm,ro,syma,zarot2
3524 : logical :: antiferro,has_qphase,my_atmtab_allocated,noncoll
3525 : logical :: paral_atom,paral_atom_unsym,use_afm,use_res
3526 : character(len=8) :: pertstrg,wrt_mode
3527 : character(len=500) :: msg
3528 : !arrays
3529 : integer,parameter :: alpha(6)=(/1,2,3,3,3,2/),beta(6)=(/1,2,3,2,1,1/)
3530 : integer :: nsym_used(2)
3531 22583 : integer, pointer :: indlmn(:,:)
3532 22583 : integer,pointer :: my_atmtab(:)
3533 22583 : integer,allocatable :: symrec_det(:)
3534 : real(dp) :: fact(2),factsym(2),phase(2),rhoijc(2),rotmag(2,3,2),rotrho(2,2,2)
3535 : real(dp) :: summag(2,3,2),sumrho(2,2,2),sum1(2),work1(2,3,3),xsym(3)
3536 22583 : real(dp),allocatable :: rotgr(:,:,:,:),rotmaggr(:,:,:,:),sumgr(:,:,:),summaggr(:,:,:,:)
3537 22583 : real(dp),allocatable :: symrec_cart(:,:,:)
3538 22583 : real(dp),pointer :: grad(:,:,:)
3539 22583 : type(coeff3_type),target,allocatable :: tmp_grhoij(:)
3540 22583 : type(pawrhoij_type),pointer :: pawrhoij_unsym_all(:)
3541 :
3542 : ! *********************************************************************
3543 :
3544 : !Sizes of pawrhoij datastructures
3545 22583 : nrhoij=size(pawrhoij)
3546 22583 : nrhoij_unsym=size(pawrhoij_unsym)
3547 :
3548 : !Set up parallelism over atoms
3549 22583 : paral_atom=(present(comm_atom).and.(nrhoij/=natom))
3550 22427 : paral_atom_unsym=(present(comm_atom).and.(nrhoij_unsym/=natom))
3551 22583 : nullify(my_atmtab);if (present(mpi_atmtab)) my_atmtab => mpi_atmtab
3552 22583 : my_comm_atom=xmpi_comm_self;if (present(comm_atom)) my_comm_atom=comm_atom
3553 22583 : call get_my_atmtab(my_comm_atom,my_atmtab,my_atmtab_allocated,paral_atom,natom)
3554 :
3555 : !Test: consistency between choice and ngrhoij
3556 22583 : ngrhoij=0
3557 22583 : if (nrhoij>0) then
3558 20805 : ngrhoij=pawrhoij(1)%ngrhoij
3559 20805 : if(choice==2) ngrhoij=min(3,pawrhoij(1)%ngrhoij)
3560 20805 : if(choice==3.or.choice==4) ngrhoij=min(6,pawrhoij(1)%ngrhoij)
3561 20805 : if(choice==23.or.choice==24) ngrhoij=min(9,pawrhoij(1)%ngrhoij)
3562 : if ((choice==1.and.ngrhoij/=0).or.(choice==2.and.ngrhoij/=3).or. &
3563 : & (choice==3.and.ngrhoij/=6).or.(choice==23.and.ngrhoij/=9).or. &
3564 20805 : & (choice==4.and.ngrhoij/=6).or.(choice==24.and.ngrhoij/=9) ) then
3565 0 : msg='Inconsistency between variables choice and ngrhoij !'
3566 0 : LIBPAW_BUG(msg)
3567 : end if
3568 : end if
3569 :
3570 : !Antiferro case ?
3571 20805 : antiferro=.false.;if (nrhoij>0) antiferro=(pawrhoij(1)%nspden==2.and.pawrhoij(1)%nsppol==1)
3572 : !Non-collinear case
3573 20805 : noncoll=.false.;if (nrhoij>0) noncoll=(pawrhoij(1)%nspden==4)
3574 : !Do we use antiferro symmetries ?
3575 22583 : use_afm=((antiferro).or.(noncoll.and.afm_noncoll))
3576 : !Do we impose zero magnetization?
3577 22583 : use_zeromag_=.false. ; if (present(use_zeromag)) use_zeromag_=use_zeromag
3578 :
3579 : ! Does not symmetrize imaginary part for GS calculations
3580 22583 : cplex_eff=1
3581 22583 : if (nrhoij>0.and.(ipert>0.or.antiferro.or.noncoll)) cplex_eff=pawrhoij(1)%cplex_rhoij
3582 :
3583 : !Do we have a phase due to q-vector?
3584 20805 : has_qphase=.false.
3585 : if (nrhoij>0) then
3586 20805 : has_qphase=(pawrhoij(1)%qphase==2)
3587 20805 : if (present(qphon)) then
3588 80779 : if (any(abs(qphon(1:3))>tol8).and.(.not.has_qphase)) then
3589 0 : msg='Should have qphase=2 for a non-zero q!'
3590 0 : LIBPAW_BUG(msg)
3591 : end if
3592 : end if
3593 : end if
3594 :
3595 : !Printing of unsymetrized Rhoij
3596 20805 : if (nrhoij>0.and.optrhoij==1.and.pawprtvol/=0) then
3597 4685 : wrt_mode='COLL';if (paral_atom) wrt_mode='PERS'
3598 4685 : pertstrg="RHOIJ";if (ipert>0) pertstrg="RHOIJ(1)"
3599 4685 : natinc=1;if(nrhoij>1.and.pawprtvol>=0) natinc=nrhoij-1
3600 4685 : write(msg, '(7a)') ch10," PAW TEST:",ch10,&
3601 9370 : & ' ========= Values of ',trim(pertstrg),' before symetrization =========',ch10
3602 4685 : call wrtout(std_out,msg,wrt_mode)
3603 13393 : do iatm=1,nrhoij,natinc
3604 8708 : iatom=iatm; if (paral_atom) iatom=my_atmtab(iatm)
3605 8708 : if (nrhoij==1.and.ipert>0.and.ipert<=natom) iatom=ipert
3606 : call pawrhoij_print_rhoij(pawrhoij_unsym(iatm)%rhoij_,pawrhoij_unsym(iatm)%cplex_rhoij,&
3607 : & pawrhoij_unsym(iatm)%qphase,iatom,natom,&
3608 13393 : & unit=std_out,opt_prtvol=pawprtvol,mode_paral=wrt_mode)
3609 : end do
3610 4685 : call wrtout(std_out,"",wrt_mode)
3611 : end if
3612 :
3613 : !Symetrization occurs only when nsym>1
3614 22583 : if (nsym>1) then
3615 :
3616 : ! Symetrization of gradients not compatible with nspden=4
3617 13227 : if (nrhoij>0) then
3618 12723 : if (choice>2.and.pawrhoij(1)%nspden==4) then
3619 0 : msg='For the time being, choice>2 is not compatible with nspden=4 !'
3620 0 : LIBPAW_BUG(msg)
3621 : end if
3622 : end if
3623 :
3624 : ! Symetry matrixes must be in memory
3625 13227 : if (pawang%nsym==0) then
3626 0 : msg='pawang%zarot must be allocated !'
3627 0 : LIBPAW_BUG(msg)
3628 : end if
3629 :
3630 13227 : if (has_qphase.and.choice>1) then
3631 0 : msg='choice>1 not compatible with q-phase !'
3632 0 : LIBPAW_BUG(msg)
3633 : end if
3634 :
3635 : ! Several inits/allocations
3636 13227 : if (noncoll) then
3637 3603 : LIBPAW_ALLOCATE(symrec_cart,(3,3,nsym))
3638 3603 : LIBPAW_ALLOCATE(symrec_det,(nsym))
3639 9149 : do irot=1,nsym
3640 7948 : symrec_cart(:,:,irot)=symrhoij_symcart(gprimd,rprimd,symrec(:,:,irot))
3641 : ! compute the sign of the determinant of the symmetries
3642 : ! to be able to apply only the proper part of the symmetries to the magn. components
3643 : ! (magnetization == pseudo-vector)
3644 : det = symrec_cart(1,1,irot)*symrec_cart(2,2,irot)*symrec_cart(3,3,irot)+&
3645 : & symrec_cart(2,1,irot)*symrec_cart(3,2,irot)*symrec_cart(1,3,irot)+&
3646 : & symrec_cart(1,2,irot)*symrec_cart(2,3,irot)*symrec_cart(3,1,irot) - &
3647 : & (symrec_cart(3,1,irot)*symrec_cart(2,2,irot)*symrec_cart(1,3,irot)+&
3648 : & symrec_cart(2,1,irot)*symrec_cart(1,2,irot)*symrec_cart(3,3,irot)+&
3649 7948 : & symrec_cart(3,2,irot)*symrec_cart(2,3,irot)*symrec_cart(1,1,irot))
3650 9149 : symrec_det(irot) = nint(det) ! should return 1 or -1
3651 : end do
3652 : end if
3653 13227 : ishift2=0;ishift3=0;ishift4=0
3654 13227 : if (choice>1) then
3655 31 : iafm=merge(2,1,antiferro)
3656 31 : qphase=merge(2,1,has_qphase)
3657 155 : LIBPAW_ALLOCATE(sumgr,(cplex_eff,ngrhoij,qphase))
3658 186 : LIBPAW_ALLOCATE(rotgr,(cplex_eff,ngrhoij,iafm,qphase))
3659 31 : if (noncoll) then
3660 0 : LIBPAW_ALLOCATE(summaggr,(cplex_eff,ngrhoij,3,qphase))
3661 0 : LIBPAW_ALLOCATE(rotmaggr,(cplex_eff,ngrhoij,3,qphase))
3662 : end if
3663 31 : if (choice==23) ishift2=6
3664 31 : if (choice==24) ishift4=3
3665 31 : if (.not.paral_atom_unsym) then
3666 : ! Have to make a temporary copy of grhoij
3667 181 : LIBPAW_DATATYPE_ALLOCATE(tmp_grhoij,(nrhoij))
3668 119 : do iatm=1,nrhoij
3669 : sz1=pawrhoij_unsym(iatm)%cplex_rhoij*pawrhoij_unsym(iatm)%qphase &
3670 88 : & *pawrhoij_unsym(iatm)%lmn2_size
3671 88 : sz2=pawrhoij_unsym(iatm)%nspden
3672 440 : LIBPAW_ALLOCATE(tmp_grhoij(iatm)%value,(ngrhoij,sz1,sz2))
3673 : tmp_grhoij(iatm)%value(1:ngrhoij,1:sz1,1:sz2)= &
3674 17723 : & pawrhoij_unsym(iatm)%grhoij(1:ngrhoij,1:sz1,1:sz2)
3675 : end do
3676 : end if
3677 : end if
3678 :
3679 : ! In case of paw_rhoij_unsym distributed over atomic sites, gather it
3680 13227 : if (paral_atom_unsym) then
3681 0 : LIBPAW_DATATYPE_ALLOCATE(pawrhoij_unsym_all,(natom))
3682 0 : call pawrhoij_nullify(pawrhoij_unsym_all)
3683 : call pawrhoij_gather(pawrhoij_unsym,pawrhoij_unsym_all,-1,my_comm_atom,&
3684 : & with_lmnmix=.false.,with_rhoijp=.false.,&
3685 0 : & with_rhoijres=.false.,with_grhoij=(choice>1))
3686 0 : nrhoij1=natom
3687 : else
3688 13227 : pawrhoij_unsym_all=>pawrhoij_unsym
3689 13227 : nrhoij1=nrhoij_unsym
3690 : end if
3691 :
3692 :
3693 : ! Loops over atoms and spin components
3694 : ! ------------------------------------
3695 43619 : do iatm=1,nrhoij
3696 30392 : iatom=iatm;if (paral_atom) iatom=my_atmtab(iatm)
3697 30392 : if (nrhoij==1.and.ipert>0.and.ipert<=natom.and.(.not.paral_atom)) iatom=ipert
3698 30392 : itypat=typat(iatom)
3699 30392 : lmn_size=pawrhoij(iatm)%lmn_size
3700 30392 : lmn2_size=pawrhoij(iatm)%lmn2_size
3701 30392 : qphase=pawrhoij(iatm)%qphase
3702 30392 : cplex_rhoij=pawrhoij(iatm)%cplex_rhoij
3703 30392 : cplex_eff=1;if (ipert>0.or.antiferro.or.noncoll) cplex_eff=cplex_rhoij
3704 30392 : use_res=(pawrhoij(iatm)%use_rhoijres>0)
3705 30392 : indlmn => pawtab(itypat)%indlmn
3706 :
3707 30392 : nselect=0
3708 64519 : do ispden=1,pawrhoij(iatm)%nsppol
3709 34127 : jspden=min(3-ispden,pawrhoij(iatm)%nsppol)
3710 :
3711 : ! Store old -rhoij in residual
3712 34127 : if (optrhoij==1.and.use_res) then
3713 2473445 : pawrhoij(iatm)%rhoijres(:,ispden)=zero
3714 1371100 : if (noncoll) pawrhoij(iatm)%rhoijres(:,2:4)=zero
3715 102342 : if (antiferro) pawrhoij(iatm)%rhoijres(:,2)=zero
3716 52994 : do iq=1,qphase
3717 26578 : iq0=(iq-1)*cplex_rhoij*lmn2_size
3718 26578 : if (cplex_rhoij==1) then
3719 964122 : do irhoij=1,pawrhoij(iatm)%nrhoijsel
3720 938796 : klmn1=iq0+pawrhoij(iatm)%rhoijselect(irhoij);jrhoij=iq0+irhoij
3721 964122 : pawrhoij(iatm)%rhoijres(klmn1,ispden)=-pawrhoij(iatm)%rhoijp(jrhoij,ispden)
3722 : end do
3723 : else
3724 167537 : do irhoij=1,pawrhoij(iatm)%nrhoijsel
3725 166285 : klmn1=iq0+2*pawrhoij(iatm)%rhoijselect(irhoij);jrhoij=iq0+2*irhoij
3726 500107 : pawrhoij(iatm)%rhoijres(klmn1-1:klmn1,ispden)=-pawrhoij(iatm)%rhoijp(jrhoij-1:jrhoij,ispden)
3727 : end do
3728 : end if
3729 26578 : if (noncoll) then
3730 1682 : if (cplex_rhoij==1) then
3731 1968 : do mu=2,4
3732 214152 : do irhoij=1,pawrhoij(iatm)%nrhoijsel
3733 212184 : klmn1=iq0+pawrhoij(iatm)%rhoijselect(irhoij);jrhoij=iq0+irhoij
3734 213660 : pawrhoij(iatm)%rhoijres(klmn1,mu)=-pawrhoij(iatm)%rhoijp(jrhoij,mu)
3735 : end do
3736 : end do
3737 : else
3738 4760 : do mu=2,4
3739 497942 : do irhoij=1,pawrhoij(iatm)%nrhoijsel
3740 493182 : klmn1=iq0+2*pawrhoij(iatm)%rhoijselect(irhoij);jrhoij=iq0+2*irhoij
3741 1483116 : pawrhoij(iatm)%rhoijres(klmn1-1:klmn1,mu)=-pawrhoij(iatm)%rhoijp(jrhoij-1:jrhoij,mu)
3742 : end do
3743 : end do
3744 : end if
3745 : end if
3746 52994 : if (antiferro) then
3747 656 : if (cplex_rhoij==1) then
3748 48512 : do irhoij=1,pawrhoij(iatm)%nrhoijsel
3749 47856 : klmn1=iq0+pawrhoij(iatm)%rhoijselect(irhoij);jrhoij=iq0+irhoij
3750 48512 : pawrhoij(iatm)%rhoijres(klmn1,2)=-pawrhoij(iatm)%rhoijp(jrhoij,2)
3751 : end do
3752 : else
3753 0 : do irhoij=1,pawrhoij(iatm)%nrhoijsel
3754 0 : klmn1=iq0+2*pawrhoij(iatm)%rhoijselect(irhoij);jrhoij=iq0+2*irhoij
3755 0 : pawrhoij(iatm)%rhoijres(klmn1-1:klmn1,2)=-pawrhoij(iatm)%rhoijp(jrhoij-1:jrhoij,2)
3756 : end do
3757 : end if
3758 : end if
3759 : end do
3760 : end if
3761 :
3762 :
3763 : ! Loops over (il,im) and (jl,jm)
3764 : ! ------------------------------
3765 : jl0=-1;jln0=-1;indexj=1
3766 431941 : do jlmn=1,lmn_size
3767 367422 : jl=indlmn(1,jlmn)
3768 367422 : jlpm=1+jl+indlmn(2,jlmn)
3769 367422 : jln=indlmn(5,jlmn)
3770 367422 : if (jln/=jln0) indexj=indexj+2*jl0+1
3771 367422 : j0lmn=jlmn*(jlmn-1)/2
3772 367422 : il0=-1;iln0=-1;indexi=1
3773 2951998 : do ilmn=1,jlmn
3774 2584576 : il=indlmn(1,ilmn)
3775 2584576 : ilpm=1+il+indlmn(2,ilmn)
3776 2584576 : iln=indlmn(5,ilmn)
3777 2584576 : if (iln/=iln0) indexi=indexi+2*il0+1
3778 2584576 : klmn=j0lmn+ilmn
3779 2584576 : klmn1=merge(klmn,2*klmn-1,cplex_rhoij==1)
3780 :
3781 2584576 : nsym_used(:)=0
3782 2584576 : if (optrhoij==1) rotrho=zero
3783 2584576 : if (optrhoij==1.and.noncoll) rotmag=zero
3784 2619600 : if (choice>1) rotgr=zero
3785 2584576 : if (choice>1.and.noncoll) rotmaggr=zero
3786 :
3787 :
3788 : ! Loop over symmetries
3789 : ! --------------------
3790 88094224 : do irot=1,nsym
3791 :
3792 85509648 : if ((symafm(irot)/=1).and.(.not.use_afm)) cycle
3793 83296272 : kspden=ispden;if (symafm(irot)==-1) kspden=jspden
3794 83296272 : iafm=1;if ((antiferro).and.(symafm(irot)==-1)) iafm=2
3795 83296272 : factafm=dble(symafm(irot))
3796 :
3797 83296272 : nsym_used(iafm)=nsym_used(iafm)+1
3798 83296272 : at_indx=min(indsym(4,irot,iatom),nrhoij1)
3799 :
3800 83296272 : if (has_qphase) then
3801 : arg=two_pi*(qphon(1)*indsym(1,irot,iatom)+qphon(2)*indsym(2,irot,iatom) &
3802 24192 : & +qphon(3)*indsym(3,irot,iatom))
3803 24192 : phase(1)=cos(arg);phase(2)=sin(arg)
3804 : end if
3805 :
3806 83296272 : if (optrhoij==1) sumrho=zero
3807 83296272 : if (optrhoij==1.and.noncoll) summag=zero
3808 84099984 : if (choice>1) sumgr=zero
3809 83296272 : if (choice>1.and.noncoll) summaggr=zero
3810 :
3811 83296272 : if (choice>1) then
3812 114816 : if (paral_atom_unsym) then
3813 0 : grad => pawrhoij_unsym_all(at_indx)%grhoij
3814 : else
3815 114816 : grad => tmp_grhoij(at_indx)%value
3816 : end if
3817 : end if
3818 :
3819 :
3820 : ! Accumulate values over (mi,mj)
3821 : ! ------------------------------
3822 439946100 : do mj=1,2*jl+1
3823 356649828 : indexjj=indexj+mj;indexjj0=indexjj*(indexjj-1)/2
3824 1627927924 : do mi=1,2*il+1
3825 3563945472 : factsym(:)=one
3826 1187981824 : indexii=indexi+mi
3827 1187981824 : if (indexii<=indexjj) then
3828 1039179472 : indexk=indexjj0+indexii
3829 1039179472 : factsym(2)=one
3830 : else
3831 148802352 : indexk=indexii*(indexii-1)/2+indexjj
3832 148802352 : factsym(2)=-one
3833 : end if
3834 1187981824 : indexkc=cplex_rhoij*(indexk-1)
3835 1187981824 : indexkc_q=indexkc+cplex_rhoij*lmn2_size
3836 :
3837 : ! Be careful: use here R_rel^-1 in term of spherical harmonics
3838 : ! which is tR_rec in term of spherical harmonics
3839 : ! so, use transpose[zarot]
3840 1187981824 : zarot2=pawang%zarot(mi,ilpm,il+1,irot)*pawang%zarot(mj,jlpm,jl+1,irot)
3841 : ! zarot2=pawang%zarot(ilpm,mi,il+1,irot)*pawang%zarot(jlpm,mj,jl+1,irot)
3842 :
3843 : ! Rotate rhoij
3844 1187981824 : if (optrhoij==1) then
3845 1186721440 : fact(1)=factsym(1);fact(2)=factsym(2)*factafm !????? What? MT
3846 : sumrho(1:cplex_eff,iafm,1)=sumrho(1:cplex_eff,iafm,1) &
3847 : & +fact(1:cplex_eff)*zarot2 &
3848 2402669180 : & *pawrhoij_unsym_all(at_indx)%rhoij_(indexkc+1:indexkc+cplex_eff,kspden)
3849 1186721440 : if (qphase==2) &
3850 : & sumrho(1:cplex_eff,iafm,2)=sumrho(1:cplex_eff,iafm,2) &
3851 : & +fact(1:cplex_eff)*zarot2 &
3852 452352 : & *pawrhoij_unsym_all(at_indx)%rhoij_(indexkc_q+1:indexkc_q+cplex_eff,kspden)
3853 : end if
3854 :
3855 : ! Rotate rhoij magnetization
3856 1187981824 : if (optrhoij==1.and.noncoll) then
3857 34442712 : fact(1)=factsym(1)*factafm;fact(2)=factsym(2)
3858 137770848 : do mu=1,3
3859 : summag(1:cplex_eff,mu,1)=summag(1:cplex_eff,mu,1) &
3860 : & +fact(1:cplex_eff)*zarot2 &
3861 328777884 : & *pawrhoij_unsym_all(at_indx)%rhoij_(indexkc+1:indexkc+cplex_eff,1+mu)
3862 : end do
3863 34442712 : if (qphase==2) then
3864 0 : do mu=1,3
3865 : summag(1:cplex_eff,mu,2)=summag(1:cplex_eff,mu,2) &
3866 : & +fact(1:cplex_eff)*zarot2 &
3867 0 : & *pawrhoij_unsym_all(at_indx)%rhoij_(indexkc_q+1:indexkc_q+cplex_eff,1+mu)
3868 : end do
3869 : end if
3870 : end if
3871 :
3872 : ! Rotate gradients of rhoij
3873 1187981824 : if (choice>1) then
3874 1260384 : fact(1)=factsym(1);fact(2)=factsym(2)*factafm !????? What? MT
3875 5041536 : do mu=1,ngrhoij
3876 : sumgr(1:cplex_eff,mu,1)=sumgr(1:cplex_eff,mu,1) &
3877 8822688 : & +fact(1:cplex_eff)*zarot2*grad(mu,indexkc+1:indexkc+cplex_eff,kspden)
3878 : end do
3879 1260384 : if (qphase==2) then
3880 0 : do mu=1,ngrhoij
3881 : sumgr(1:cplex_eff,mu,2)=sumgr(1:cplex_eff,mu,2) &
3882 0 : & +fact(1:cplex_eff)*zarot2*grad(mu,indexkc_q+1:indexkc_q+cplex_eff,kspden)
3883 : end do
3884 : end if
3885 : end if
3886 :
3887 : ! Rotate gradients of rhoij magnetization
3888 1544631652 : if (choice>1.and.noncoll) then
3889 0 : fact(1)=factsym(1)*factafm;fact(2)=factsym(2)
3890 0 : do mu=1,3
3891 0 : do nu=1,ngrhoij
3892 : summaggr(1:cplex_eff,nu,mu,1)=summaggr(1:cplex_eff,nu,mu,1) &
3893 0 : & +fact(1:cplex_eff)*zarot2*grad(nu,indexkc+1:indexkc+cplex_eff,1+mu)
3894 : end do
3895 : end do
3896 0 : if (qphase==2) then
3897 0 : do mu=1,3
3898 0 : do nu=1,ngrhoij
3899 : summaggr(1:cplex_eff,nu,mu,2)=summaggr(1:cplex_eff,nu,mu,2) &
3900 0 : & +fact(1:cplex_eff)*zarot2*grad(nu,indexkc_q+1:indexkc_q+cplex_eff,1+mu)
3901 : end do
3902 : end do
3903 : end if
3904 : end if
3905 :
3906 : end do ! mi
3907 : end do ! mj
3908 :
3909 : ! Apply phase for phonons
3910 83296272 : if (has_qphase) then
3911 : !Remember, RHOij is stored as follows:
3912 : ! RHOij= [rhoij(2klmn-1)+i.rhoij(2klmn)]
3913 : ! +i.[rhoij(2lnm2_size+2klmn-1)+i.rhoij(2lmn2_size+2klmn)]
3914 24192 : if (optrhoij==1) then
3915 48384 : do iplex=1,cplex_rhoij
3916 24192 : rhoijc(1)=sumrho(iplex,iafm,1)
3917 24192 : rhoijc(2)=sumrho(iplex,iafm,2)
3918 24192 : sumrho(iplex,iafm,1)=phase(1)*rhoijc(1)-phase(2)*rhoijc(2)
3919 48384 : sumrho(iplex,iafm,2)=phase(1)*rhoijc(2)+phase(2)*rhoijc(1)
3920 : end do
3921 24192 : if (noncoll) then
3922 0 : do iplex=1,cplex_rhoij
3923 0 : do mu=1,3
3924 0 : rhoijc(1)=summag(iplex,mu,1)
3925 0 : rhoijc(2)=summag(iplex,mu,2)
3926 0 : summag(iplex,mu,1)=phase(1)*rhoijc(1)-phase(2)*rhoijc(2)
3927 0 : summag(iplex,mu,2)=phase(1)*rhoijc(2)+phase(2)*rhoijc(1)
3928 : end do
3929 : end do
3930 : end if
3931 : end if
3932 24192 : if (choice>1) then
3933 0 : do iplex=1,cplex_rhoij
3934 0 : do mu=1,3
3935 0 : rhoijc(1)=sumgr(iplex,mu,1)
3936 0 : rhoijc(2)=sumgr(iplex,mu,2)
3937 0 : sumgr(iplex,mu,1)=phase(1)*rhoijc(1)-phase(2)*rhoijc(2)
3938 0 : sumgr(iplex,mu,2)=phase(1)*rhoijc(2)+phase(2)*rhoijc(1)
3939 : end do
3940 : end do
3941 0 : if (noncoll) then
3942 0 : do iplex=1,cplex_rhoij
3943 0 : do mu=1,3
3944 0 : do nu=1,ngrhoij
3945 0 : rhoijc(1)=summaggr(iplex,nu,mu,1)
3946 0 : rhoijc(2)=summaggr(iplex,nu,mu,2)
3947 0 : summaggr(iplex,nu,mu,1)=phase(1)*rhoijc(1)-phase(2)*rhoijc(2)
3948 0 : summaggr(iplex,nu,mu,2)=phase(1)*rhoijc(2)+phase(2)*rhoijc(1)
3949 : end do
3950 : end do
3951 : end do
3952 : end if
3953 : end if
3954 : end if
3955 :
3956 : ! Add contribution of this rotation
3957 83296272 : if (optrhoij==1) then
3958 166387104 : do iq=1,qphase
3959 : rotrho(1:cplex_eff,iafm,iq)=rotrho(1:cplex_eff,iafm,iq) &
3960 251141132 : & +sumrho(1:cplex_eff,iafm,iq)
3961 : end do
3962 : end if
3963 :
3964 :
3965 : ! Rotate vector fields in real space (forces, magnetization, etc...)
3966 : ! Should use symrel^-1 but use transpose[symrec] instead
3967 : ! ===== Rhoij magnetization ====
3968 83296272 : if (noncoll.and.optrhoij==1) then
3969 3782128 : do iq=1,qphase
3970 9455320 : do nu=1,3
3971 24583832 : do mu=1,3
3972 : rotmag(1:cplex_eff,mu,iq)=rotmag(1:cplex_eff,mu,iq) &
3973 53647764 : & +symrec_det(irot)*symrec_cart(mu,nu,irot)*summag(1:cplex_eff,nu,iq)
3974 : end do
3975 : end do
3976 : end do
3977 : end if
3978 : ! ===== Derivatives vs atomic positions ====
3979 83296272 : if (choice==2.or.choice==23.or.choice==24) then
3980 229632 : do iq=1,qphase
3981 574080 : do nu=1,3
3982 344448 : nushift=nu+ishift2
3983 1492608 : do mu=1,3
3984 1033344 : mushift=mu+ishift2
3985 : rotgr(1:cplex_eff,mushift,iafm,iq)=rotgr(1:cplex_eff,mushift,iafm,iq) &
3986 2411136 : & +dble(symrec(mu,nu,irot))*sumgr(1:cplex_eff,nushift,iq)
3987 : end do
3988 : end do
3989 : end do
3990 114816 : if (noncoll) then
3991 0 : do iq=1,qphase
3992 0 : do mub=1,3 ! Loop on magnetization components
3993 0 : do mua=1,3 ! Loop on gradients
3994 0 : mushift=mua+ishift2
3995 0 : sum1(:)=zero;xsym(1:3)=dble(symrec(mua,1:3,irot))
3996 0 : do nu=1,3
3997 0 : syma=symrec_det(irot)*symrec_cart(mub,nu,irot)
3998 : sum1(1:cplex_eff)=sum1(1:cplex_eff)+syma &
3999 : & *(summaggr(1:cplex_eff,ishift2+1,nu,iq)*xsym(1) &
4000 : & +summaggr(1:cplex_eff,ishift2+2,nu,iq)*xsym(2) &
4001 0 : & +summaggr(1:cplex_eff,ishift2+3,nu,iq)*xsym(3))
4002 : end do
4003 : rotmaggr(1:cplex_eff,mushift,mub,iq)= &
4004 0 : & rotmaggr(1:cplex_eff,mushift,mub,iq)+sum1(1:cplex_eff)
4005 : end do
4006 : end do
4007 : end do
4008 : end if
4009 : end if
4010 : ! ===== Derivatives vs strain ====
4011 83296272 : if (choice==3.or.choice==23) then
4012 0 : do iq=1,qphase
4013 0 : work1(1:cplex_eff,1,1)=sumgr(1:cplex_eff,1+ishift3,iq);work1(1:cplex_eff,2,2)=sumgr(1:cplex_eff,2+ishift3,iq)
4014 0 : work1(1:cplex_eff,3,3)=sumgr(1:cplex_eff,3+ishift3,iq);work1(1:cplex_eff,2,3)=sumgr(1:cplex_eff,4+ishift3,iq)
4015 0 : work1(1:cplex_eff,1,3)=sumgr(1:cplex_eff,5+ishift3,iq);work1(1:cplex_eff,1,2)=sumgr(1:cplex_eff,6+ishift3,iq)
4016 0 : work1(1:cplex_eff,3,1)=work1(1:cplex_eff,1,3);work1(1:cplex_eff,3,2)=work1(1:cplex_eff,2,3)
4017 0 : work1(1:cplex_eff,2,1)=work1(1:cplex_eff,1,2)
4018 0 : do mu=1,6
4019 0 : mushift=mu+ishift3
4020 0 : mua=alpha(mu);mub=beta(mu)
4021 0 : sum1(:)=zero;xsym(1:3)=dble(symrec(mub,1:3,irot))
4022 0 : do nu=1,3
4023 0 : syma=dble(symrec(mua,nu,irot))
4024 : sum1(1:cplex_eff)=sum1(1:cplex_eff) &
4025 : & +syma*(work1(1:cplex_eff,nu,1)*xsym(1) &
4026 : & +work1(1:cplex_eff,nu,2)*xsym(2) &
4027 0 : & +work1(1:cplex_eff,nu,3)*xsym(3))
4028 : end do
4029 : rotgr(1:cplex_eff,mushift,iafm,iq)= &
4030 0 : & rotgr(1:cplex_eff,mushift,iafm,iq)+sum1(1:cplex_eff)
4031 : end do
4032 : end do
4033 : end if
4034 : ! ===== Second derivatives vs atomic positions ====
4035 85880848 : if (choice==4.or.choice==24) then
4036 0 : do iq=1,qphase
4037 0 : work1(1:cplex_eff,1,1)=sumgr(1:cplex_eff,1+ishift4,iq);work1(1:cplex_eff,2,2)=sumgr(1:cplex_eff,2+ishift4,iq)
4038 0 : work1(1:cplex_eff,3,3)=sumgr(1:cplex_eff,3+ishift4,iq);work1(1:cplex_eff,2,3)=sumgr(1:cplex_eff,4+ishift4,iq)
4039 0 : work1(1:cplex_eff,1,3)=sumgr(1:cplex_eff,5+ishift4,iq);work1(1:cplex_eff,1,2)=sumgr(1:cplex_eff,6+ishift4,iq)
4040 0 : work1(1:cplex_eff,3,1)=work1(1:cplex_eff,1,3);work1(1:cplex_eff,3,2)=work1(1:cplex_eff,2,3)
4041 0 : work1(1:cplex_eff,2,1)=work1(1:cplex_eff,1,2)
4042 0 : do mu=1,6
4043 0 : mushift=mu+ishift4
4044 0 : mua=alpha(mu);mub=beta(mu)
4045 0 : sum1(:)=zero
4046 0 : xsym(1:3)=dble(symrec(mub,1:3,irot))
4047 0 : do nu=1,3
4048 0 : syma=dble(symrec(mua,nu,irot))
4049 : sum1(1:cplex_eff)=sum1(1:cplex_eff) &
4050 : & +syma*(work1(1:cplex_eff,nu,1)*xsym(1) &
4051 : & +work1(1:cplex_eff,nu,2)*xsym(2) &
4052 0 : & +work1(1:cplex_eff,nu,3)*xsym(3))
4053 : end do
4054 : rotgr(1:cplex_eff,mushift,iafm,iq)= &
4055 0 : & rotgr(1:cplex_eff,mushift,iafm,iq)+sum1(1:cplex_eff)
4056 : end do
4057 : end do
4058 : end if
4059 :
4060 : end do ! End loop over symmetries
4061 :
4062 :
4063 : ! Store average result (over symmetries)
4064 : ! --------------------------------------
4065 :
4066 : ! Rhoij
4067 2584576 : if (optrhoij==1) then
4068 5172492 : do iq=1,qphase
4069 2592294 : klmn1q=klmn1+(iq-1)*lmn2_size*cplex_rhoij
4070 2592294 : pawrhoij(iatm)%rhoijp(klmn1q,ispden)=rotrho(1,1,iq)/nsym_used(1)
4071 5172492 : if (cplex_rhoij==2) then
4072 183439 : if (cplex_eff==1) ro=pawrhoij_unsym_all(iatom)%rhoij_(klmn1q+1,ispden)
4073 183439 : if (cplex_eff==2) ro=rotrho(2,1,iq)/nsym_used(1)
4074 183439 : pawrhoij(iatm)%rhoijp(klmn1q+1,ispden)=ro
4075 : end if
4076 : end do
4077 : end if
4078 :
4079 : ! Rhoij magnetization
4080 2584576 : if (noncoll.and.optrhoij==1) then
4081 1061356 : do mu=2,4
4082 1857373 : do iq=1,qphase
4083 796017 : klmn1q=klmn1+(iq-1)*lmn2_size*cplex_rhoij
4084 796017 : if (use_zeromag_) then
4085 0 : pawrhoij(iatm)%rhoijp(klmn1q,mu)=zero
4086 : else
4087 796017 : pawrhoij(iatm)%rhoijp(klmn1q,mu)=rotmag(1,mu-1,iq)/nsym_used(1)
4088 : end if
4089 1592034 : if (cplex_rhoij==2) then
4090 543621 : if (cplex_eff==1) ro=pawrhoij_unsym_all(iatom)%rhoij_(klmn1q+1,mu)
4091 543621 : if (cplex_eff==2) ro=rotmag(2,mu-1,iq)/nsym_used(1)
4092 543621 : pawrhoij(iatm)%rhoijp(klmn1q+1,mu)=ro
4093 : end if
4094 : end do
4095 : end do
4096 : end if
4097 :
4098 : ! Rhoij^down when antiferro
4099 2584576 : if (antiferro.and.optrhoij==1) then
4100 75926 : if (nsym_used(2)>0) then
4101 151852 : do iq=1,qphase
4102 75926 : klmn1q=klmn1+(iq-1)*lmn2_size*cplex_rhoij
4103 75926 : pawrhoij(iatm)%rhoijp(klmn1q,2)=rotrho(1,2,iq)/nsym_used(2)
4104 151852 : if (cplex_rhoij==2) then
4105 0 : if (cplex_eff==1) ro=pawrhoij_unsym_all(iatom)%rhoij_(klmn1q+1,2)
4106 0 : if (cplex_eff==2) ro=rotrho(2,2,iq)/nsym_used(2)
4107 0 : pawrhoij(iatm)%rhoijp(klmn1q+1,2)=ro
4108 : end if
4109 : end do
4110 : end if
4111 : end if
4112 :
4113 : ! Gradients of rhoij
4114 2584576 : if (choice>1) then
4115 8756 : do iq=1,qphase
4116 4378 : klmn1q=klmn1+(iq-1)*lmn2_size*cplex_rhoij
4117 13134 : do iplex=1,cplex_eff
4118 17512 : do mu=1,ngrhoij
4119 17512 : pawrhoij(iatm)%grhoij(mu,klmn1q,ispden)=rotgr(iplex,mu,1,iq)/nsym_used(1)
4120 : end do
4121 4378 : if (noncoll) then
4122 0 : if (use_zeromag_.and.iplex==1) then
4123 0 : pawrhoij(iatm)%grhoij(mu,klmn1q,2:4)=zero
4124 : else
4125 0 : do nu=1,3
4126 0 : pawrhoij(iatm)%grhoij(mu,klmn1q,1+nu)=rotmaggr(iplex,mu,nu,iq)/nsym_used(1)
4127 : end do
4128 : end if
4129 : end if
4130 4378 : if (antiferro.and.nsym_used(2)>0) then
4131 0 : do mu=1,ngrhoij
4132 0 : pawrhoij(iatm)%grhoij(mu,klmn1q,ispden)=rotgr(iplex,mu,2,iq)/nsym_used(2)
4133 : end do
4134 : end if
4135 8756 : klmn1q=klmn1q+1
4136 : end do
4137 : !if cplex_eff<cplex_rhoij, imaginary part of grhoij is unchanged
4138 : end do
4139 : end if
4140 :
4141 :
4142 2951998 : il0=il;iln0=iln ! End loops over (il,im) and (jl,jm)
4143 : end do
4144 401549 : jl0=jl;jln0=jln
4145 : end do
4146 :
4147 : end do ! End loop over ispden
4148 :
4149 : ! Select non-zero elements of rhoij
4150 30392 : if (optrhoij==1) then
4151 : call pawrhoij_filter(pawrhoij(iatm)%rhoijp,pawrhoij(iatm)%rhoijselect,&
4152 : & pawrhoij(iatm)%nrhoijsel,cplex_rhoij,qphase,lmn2_size,&
4153 30304 : & pawrhoij(iatm)%nspden)
4154 : end if
4155 :
4156 : ! Add new rhoij to rhoij residual
4157 43619 : if (optrhoij==1.and.use_res) then
4158 54811 : do ispden=1,pawrhoij(iatm)%nspden
4159 87091 : do iq=1,qphase
4160 32280 : iq0=(iq-1)*lmn2_size*cplex_rhoij
4161 64398 : if (cplex_rhoij==1) then
4162 1287208 : do irhoij=1,pawrhoij(iatm)%nrhoijsel
4163 1259750 : klmn1=iq0+pawrhoij(iatm)%rhoijselect(irhoij) ; jrhoij=iq0+irhoij
4164 : pawrhoij(iatm)%rhoijres(klmn1,ispden)= &
4165 : & pawrhoij(iatm)%rhoijres(klmn1,ispden) &
4166 1287208 : & +pawrhoij(iatm)%rhoijp(jrhoij,ispden)
4167 : end do
4168 : else
4169 672273 : do irhoij=1,pawrhoij(iatm)%nrhoijsel
4170 667451 : klmn1=iq0+2*pawrhoij(iatm)%rhoijselect(irhoij) ; jrhoij=iq0+2*irhoij
4171 : pawrhoij(iatm)%rhoijres(klmn1-1:klmn1,ispden)= &
4172 : & pawrhoij(iatm)%rhoijres(klmn1-1:klmn1,ispden) &
4173 2007175 : & +pawrhoij(iatm)%rhoijp(jrhoij-1:jrhoij,ispden)
4174 : end do
4175 : end if
4176 : end do
4177 : end do
4178 : end if
4179 :
4180 : end do ! End loop over atoms
4181 :
4182 13227 : if (noncoll) then
4183 1201 : LIBPAW_DEALLOCATE(symrec_cart)
4184 1201 : LIBPAW_DEALLOCATE(symrec_det)
4185 : end if
4186 13227 : if (choice>1) then
4187 31 : if (.not.paral_atom_unsym) then
4188 119 : do iatm=1,nrhoij
4189 119 : LIBPAW_DEALLOCATE(tmp_grhoij(iatm)%value)
4190 : end do
4191 119 : LIBPAW_DATATYPE_DEALLOCATE(tmp_grhoij)
4192 : end if
4193 31 : LIBPAW_DEALLOCATE(sumgr)
4194 31 : LIBPAW_DEALLOCATE(rotgr)
4195 31 : if (noncoll) then
4196 0 : LIBPAW_DEALLOCATE(summaggr)
4197 0 : LIBPAW_DEALLOCATE(rotmaggr)
4198 : end if
4199 : end if
4200 13227 : if(paral_atom_unsym) then
4201 0 : call pawrhoij_free(pawrhoij_unsym_all)
4202 0 : LIBPAW_DATATYPE_DEALLOCATE(pawrhoij_unsym_all)
4203 : end if
4204 :
4205 :
4206 : else ! nsym>1
4207 :
4208 : ! *********************************************************************
4209 : ! If nsym==1, only copy rhoij_ into rhoij
4210 : ! also has to fill rhoijselect array
4211 :
4212 9356 : if (antiferro) then
4213 0 : msg=' In the antiferromagnetic case, nsym cannot be 1'
4214 0 : LIBPAW_BUG(msg)
4215 : end if
4216 :
4217 9356 : if (optrhoij==1) then
4218 :
4219 28718 : do iatm=1,nrhoij
4220 19363 : iatom=iatm;if ((paral_atom).and.(.not.paral_atom_unsym)) iatom=my_atmtab(iatm)
4221 19363 : cplex_rhoij=pawrhoij(iatm)%cplex_rhoij
4222 19363 : qphase=pawrhoij(iatm)%qphase
4223 19363 : lmn2_size=pawrhoij(iatm)%lmn2_size
4224 19363 : use_res=(pawrhoij(iatm)%use_rhoijres>0)
4225 :
4226 : ! Store -rhoij_input in rhoij residual
4227 19363 : if (use_res) then
4228 922187 : pawrhoij(iatm)%rhoijres(:,:)=zero
4229 33974 : do iq=1,qphase
4230 17191 : iq0=(iq-1)*lmn2_size*cplex_rhoij
4231 33974 : if (cplex_rhoij==1) then
4232 32690 : do ispden=1,pawrhoij(iatm)%nspden
4233 381611 : do irhoij=1,pawrhoij(iatm)%nrhoijsel
4234 348921 : klmn=iq0+pawrhoij(iatm)%rhoijselect(irhoij);jrhoij=iq0+irhoij
4235 365828 : pawrhoij(iatm)%rhoijres(klmn,ispden)=-pawrhoij(iatm)%rhoijp(jrhoij,ispden)
4236 : end do
4237 : end do
4238 : else
4239 6115 : do ispden=1,pawrhoij(iatm)%nspden
4240 231127 : do irhoij=1,pawrhoij(iatm)%nrhoijsel
4241 225012 : klmn1=iq0+2*pawrhoij(iatm)%rhoijselect(irhoij);jrhoij=iq0+2*irhoij
4242 225012 : pawrhoij(iatm)%rhoijres(klmn1-1,ispden)=-pawrhoij(iatm)%rhoijp(jrhoij-1,ispden)
4243 229719 : pawrhoij(iatm)%rhoijres(klmn1 ,ispden)=-pawrhoij(iatm)%rhoijp(jrhoij ,ispden)
4244 : end do
4245 : end do
4246 : end if
4247 : end do
4248 : end if
4249 :
4250 : ! Select non-zero elements of input rhoij
4251 : call pawrhoij_filter(pawrhoij(iatm)%rhoijp,pawrhoij(iatm)%rhoijselect,&
4252 : & pawrhoij(iatm)%nrhoijsel,cplex_rhoij,qphase,lmn2_size,&
4253 19363 : & pawrhoij(iatm)%nspden,rhoij_input=pawrhoij_unsym(iatom)%rhoij_)
4254 :
4255 : ! Add new rhoij to rhoij residual
4256 28718 : if (use_res) then
4257 37983 : do ispden=1,pawrhoij(iatm)%nspden
4258 59597 : do iq=1,qphase
4259 21614 : iq0=(iq-1)*lmn2_size*cplex_rhoij
4260 42814 : if (cplex_rhoij==1) then
4261 393309 : do irhoij=1,pawrhoij(iatm)%nrhoijsel
4262 376402 : klmn1=iq0+pawrhoij(iatm)%rhoijselect(irhoij) ; jrhoij=iq0+irhoij
4263 : pawrhoij(iatm)%rhoijres(klmn1,ispden)= &
4264 : & pawrhoij(iatm)%rhoijres(klmn1,ispden) &
4265 393309 : & +pawrhoij(iatm)%rhoijp(jrhoij,ispden)
4266 : end do
4267 : else
4268 248046 : do irhoij=1,pawrhoij(iatm)%nrhoijsel
4269 243339 : klmn1=iq0+2*pawrhoij(iatm)%rhoijselect(irhoij) ; jrhoij=iq0+2*irhoij
4270 : pawrhoij(iatm)%rhoijres(klmn1-1:klmn1,ispden)= &
4271 : & pawrhoij(iatm)%rhoijres(klmn1-1:klmn1,ispden) &
4272 734724 : & +pawrhoij(iatm)%rhoijp(jrhoij-1:jrhoij,ispden)
4273 : end do
4274 : end if
4275 : end do
4276 : end do
4277 : end if
4278 :
4279 : end do ! iatm
4280 : end if ! optrhoij
4281 :
4282 : end if
4283 :
4284 : !*********************************************************************
4285 : !Printing of symetrized Rhoij
4286 22583 : if (nrhoij>0.and.optrhoij==1.and.pawprtvol/=0) then
4287 4685 : wrt_mode='COLL';if (paral_atom) wrt_mode='PERS'
4288 4685 : pertstrg="RHOIJ";if (ipert>0) pertstrg="RHOIJ(1)"
4289 4685 : natinc=1;if(nrhoij>1.and.pawprtvol>=0) natinc=nrhoij-1
4290 4685 : write(msg, '(7a)') ch10," PAW TEST:",ch10,&
4291 9370 : & ' ========= Values of ',trim(pertstrg),' after symetrization =========',ch10
4292 4685 : call wrtout(std_out,msg,wrt_mode)
4293 13393 : do iatm=1,nrhoij,natinc
4294 8708 : iatom=iatm; if (paral_atom) iatom=my_atmtab(iatm)
4295 8708 : if (nrhoij==1.and.ipert>0.and.ipert<=natom) iatom=ipert
4296 : call pawrhoij_print_rhoij(pawrhoij(iatm)%rhoijp,pawrhoij(iatm)%cplex_rhoij,&
4297 : & pawrhoij(iatm)%qphase,iatom,natom,&
4298 : & rhoijselect=pawrhoij(iatm)%rhoijselect,unit=std_out,&
4299 13393 : & opt_prtvol=pawprtvol,mode_paral=wrt_mode)
4300 : end do
4301 4685 : call wrtout(std_out,"",wrt_mode)
4302 : end if
4303 :
4304 : !Destroy atom table used for parallelism
4305 45166 : call free_my_atmtab(my_atmtab,my_atmtab_allocated)
4306 :
4307 : !*********************************************************************
4308 : !Small function: convert a symmetry operation
4309 : !from reduced coordinates (integers) to cartesian coordinates (reals)
4310 : contains
4311 7948 : function symrhoij_symcart(aprim,bprim,symred)
4312 :
4313 : real(dp) :: symrhoij_symcart(3,3)
4314 : integer,intent(in) :: symred(3,3)
4315 : real(dp),intent(in) :: aprim(3,3),bprim(3,3)
4316 : integer :: ii,jj,kk
4317 : real(dp) :: tmp(3,3)
4318 103324 : symrhoij_symcart=zero;tmp=zero
4319 31792 : do kk=1,3
4320 103324 : do jj=1,3
4321 309972 : do ii=1,3
4322 286128 : tmp(ii,jj)=tmp(ii,jj)+bprim(ii,kk)*dble(symred(jj,kk))
4323 : end do
4324 : end do
4325 : end do
4326 31792 : do kk=1,3
4327 103324 : do jj=1,3
4328 309972 : do ii=1,3
4329 286128 : symrhoij_symcart(ii,jj)=symrhoij_symcart(ii,jj)+aprim(ii,kk)*tmp(jj,kk)
4330 : end do
4331 : end do
4332 : end do
4333 : end function symrhoij_symcart
4334 :
4335 : end subroutine pawrhoij_symrhoij
4336 : !!***
4337 :
4338 : !----------------------------------------------------------------------
4339 :
4340 : !!****f* m_pawrhoij/pawrhoij_isendreceive_getbuffer
4341 : !! NAME
4342 : !! pawrhoij_isendreceive_getbuffer
4343 : !!
4344 : !! FUNCTION
4345 : !! Fill a pawrhoij structure with the buffers received in a receive operation
4346 : !! This buffer should have been first extracted by a call to pawrhoij_isendreceive_fillbuffer
4347 : !!
4348 : !! INPUTS
4349 : !! atm_indx_recv(1:total number of atoms)= array for receive operation
4350 : !! Given an index of atom in global numbering, give its index
4351 : !! in the table of atoms treated by current processor
4352 : !! or -1 if the atoms is not treated by current processor
4353 : !! buf_int= buffer of receive integers
4354 : !! buf_dp= buffer of receive double precision numbers
4355 : !! nrhoij_send= number of sent atoms
4356 : !!
4357 : !! OUTPUT
4358 : !! pawrhoij= output datastructure filled with buffers receive in a receive operation
4359 : !!
4360 : !! SOURCE
4361 :
4362 8 : subroutine pawrhoij_isendreceive_getbuffer(pawrhoij,nrhoij_send,atm_indx_recv,buf_int,buf_dp)
4363 :
4364 : !Arguments ------------------------------------
4365 : !scalars
4366 : integer,intent(in) :: nrhoij_send
4367 : !arrays
4368 : integer,intent(in) ::atm_indx_recv(:),buf_int(:)
4369 : real(dp),intent(in) :: buf_dp(:)
4370 : type(pawrhoij_type),target,intent(inout) :: pawrhoij(:)
4371 :
4372 : !Local variables-------------------------------
4373 : !scalars
4374 : integer :: buf_dp_size,buf_int_size,cplex,ii,indx_int,indx_dp,iatom_tot,irhoij_send
4375 : integer :: isp,jj,jrhoij,lmn2_size,lmnmix,ngrhoij,nselect,nspden,qphase,rhoij_size2
4376 : integer :: use_rhoijp,use_rhoijres,use_rhoij_
4377 : character(len=500) :: msg
4378 : type(pawrhoij_type),pointer :: pawrhoij1
4379 : !arrays
4380 :
4381 : ! *********************************************************************
4382 :
4383 8 : buf_int_size=size(buf_int)
4384 8 : buf_dp_size=size(buf_dp)
4385 8 : indx_int=1;indx_dp=1
4386 :
4387 16 : do irhoij_send=1,nrhoij_send
4388 :
4389 8 : iatom_tot=buf_int(indx_int) ; indx_int=indx_int+1
4390 8 : jrhoij= atm_indx_recv(iatom_tot)
4391 8 : if (jrhoij==-1) then
4392 0 : msg="Error in pawrhoij_isendreceive_getbuffer atom not found"
4393 0 : LIBPAW_BUG(msg)
4394 : end if
4395 8 : pawrhoij1=>pawrhoij(jrhoij)
4396 :
4397 8 : cplex =buf_int(indx_int) ;indx_int=indx_int+1
4398 8 : qphase =buf_int(indx_int) ;indx_int=indx_int+1
4399 8 : lmn2_size =buf_int(indx_int) ;indx_int=indx_int+1
4400 8 : nspden =buf_int(indx_int) ;indx_int=indx_int+1
4401 8 : nselect =buf_int(indx_int) ;indx_int=indx_int+1
4402 8 : lmnmix =buf_int(indx_int) ;indx_int=indx_int+1
4403 8 : ngrhoij =buf_int(indx_int) ;indx_int=indx_int+1
4404 8 : use_rhoijp =buf_int(indx_int) ;indx_int=indx_int+1
4405 8 : use_rhoijres=buf_int(indx_int) ;indx_int=indx_int+1
4406 8 : use_rhoij_ =buf_int(indx_int) ;indx_int=indx_int+1
4407 8 : rhoij_size2 =buf_int(indx_int) ;indx_int=indx_int+1
4408 8 : pawrhoij1%itypat=buf_int(indx_int) ;indx_int=indx_int+1
4409 8 : pawrhoij1%lmn_size=buf_int(indx_int) ;indx_int=indx_int+1
4410 8 : pawrhoij1%nsppol=buf_int(indx_int) ;indx_int=indx_int+1
4411 8 : pawrhoij1%nspinor=buf_int(indx_int) ;indx_int=indx_int+1
4412 8 : pawrhoij1%cplex_rhoij=cplex
4413 8 : pawrhoij1%qphase=qphase
4414 8 : pawrhoij1%lmn2_size=lmn2_size
4415 8 : pawrhoij1%nspden=nspden
4416 8 : pawrhoij1%nrhoijsel=nselect
4417 8 : pawrhoij1%lmnmix_sz=lmnmix
4418 8 : pawrhoij1%ngrhoij=ngrhoij
4419 8 : pawrhoij1%use_rhoijp=use_rhoijp
4420 8 : pawrhoij1%use_rhoijres=use_rhoijres
4421 8 : pawrhoij1%use_rhoij_=use_rhoij_
4422 8 : if (use_rhoijp>0) then
4423 24 : LIBPAW_ALLOCATE(pawrhoij1%rhoijselect,(lmn2_size))
4424 104 : pawrhoij1%rhoijselect(1:nselect)=buf_int(indx_int:indx_int+nselect-1)
4425 200 : if (nselect < lmn2_size )pawrhoij1%rhoijselect(nselect+1:lmn2_size)=zero
4426 8 : indx_int=indx_int+nselect
4427 32 : LIBPAW_ALLOCATE(pawrhoij1%rhoijp,(cplex*qphase*lmn2_size,nspden))
4428 16 : do isp=1,nspden
4429 24 : do ii=1,qphase
4430 8 : jj=(ii-1)*cplex*lmn2_size
4431 104 : pawrhoij1%rhoijp(jj+1:jj+cplex*nselect,isp)=buf_dp(indx_dp:indx_dp+cplex*nselect-1)
4432 200 : if (nselect<lmn2_size)pawrhoij1%rhoijp(jj+cplex*nselect+1:jj+cplex*lmn2_size,isp)=zero
4433 16 : indx_dp=indx_dp+cplex*nselect
4434 : end do
4435 : end do
4436 : end if
4437 8 : if (lmnmix>0) then
4438 0 : LIBPAW_ALLOCATE(pawrhoij1%kpawmix,(lmnmix))
4439 0 : pawrhoij1%kpawmix(1:lmnmix)=buf_int(indx_int:indx_int+lmnmix-1)
4440 : indx_int=indx_int+lmnmix
4441 : end if
4442 8 : if (ngrhoij>0) then
4443 0 : LIBPAW_ALLOCATE(pawrhoij1%grhoij,(ngrhoij,cplex*qphase*lmn2_size,nspden))
4444 0 : do isp=1,nspden
4445 0 : do ii=1,cplex*qphase*lmn2_size
4446 0 : pawrhoij1%grhoij(1:ngrhoij,ii,isp)=buf_dp(indx_dp:indx_dp+ngrhoij-1)
4447 0 : indx_dp=indx_dp+ngrhoij
4448 : end do
4449 : end do
4450 : end if
4451 8 : if (use_rhoijres>0) then
4452 0 : LIBPAW_ALLOCATE(pawrhoij1%rhoijres,(cplex*qphase*lmn2_size,nspden))
4453 0 : do isp=1,nspden
4454 0 : pawrhoij1%rhoijres(1:cplex*qphase*lmn2_size,isp)=buf_dp(indx_dp:indx_dp+cplex*qphase*lmn2_size-1)
4455 0 : indx_dp=indx_dp+cplex*qphase*lmn2_size
4456 : end do
4457 : end if
4458 16 : if (use_rhoij_>0) then
4459 0 : LIBPAW_ALLOCATE(pawrhoij1%rhoij_,(cplex*qphase*lmn2_size,rhoij_size2))
4460 0 : do isp=1,rhoij_size2
4461 0 : pawrhoij1%rhoij_(1:cplex*qphase*lmn2_size,isp)=buf_dp(indx_dp:indx_dp+cplex*qphase*lmn2_size-1)
4462 0 : indx_dp=indx_dp+cplex*qphase*lmn2_size
4463 : end do
4464 : end if
4465 : end do !irhoij_send
4466 8 : if ((indx_int/=1+buf_int_size).or.(indx_dp/=1+buf_dp_size)) then
4467 0 : write(msg,'(a,i10,a,i10)') 'Wrong buffer sizes: buf_int_size=',buf_int_size,' buf_dp_size=',buf_dp_size
4468 0 : LIBPAW_BUG(msg)
4469 : end if
4470 :
4471 8 : end subroutine pawrhoij_isendreceive_getbuffer
4472 : !!***
4473 :
4474 : !----------------------------------------------------------------------
4475 :
4476 : !!****f* m_pawrhoij/pawrhoij_isendreceive_fillbuffer
4477 : !! NAME
4478 : !! pawrhoij_isendreceive_fillbuffer
4479 : !!
4480 : !! FUNCTION
4481 : !! Extract from pawrhoij and from the global index of atoms
4482 : !! the buffers to send in a sending operation
4483 : !! This function has to be coupled with a call to pawrhoij_isendreceive_getbuffer
4484 : !!
4485 : !! INPUTS
4486 : !! atm_indx_send(1:total number of atoms)= array for send operation,
4487 : !! Given an index of atom in global numbering, give its index
4488 : !! in the table of atoms treated by current processor
4489 : !! or -1 if the atoms is not treated by current processor
4490 : !! nrhoij_send= number of sent atoms
4491 : !! pawrhoij= data structure from which are extract buffer int and buffer dp
4492 : !!
4493 : !! OUTPUT
4494 : !! buf_int= buffer of integers to be sent
4495 : !! buf_int_size= size of buffer of integers
4496 : !! buf_dp= buffer of double precision numbers to be sent
4497 : !! buf_dp_size= size of buffer of double precision numbers
4498 : !!
4499 : !! SOURCE
4500 : !!
4501 8 : subroutine pawrhoij_isendreceive_fillbuffer(pawrhoij,atmtab_send, atm_indx_send,nrhoij_send,&
4502 : & buf_int,buf_int_size,buf_dp,buf_dp_size)
4503 :
4504 : !Arguments ------------------------------------
4505 : !scalars
4506 : integer,intent(out) :: buf_dp_size,buf_int_size
4507 : integer,intent(in) :: nrhoij_send
4508 : !arrays
4509 : integer,intent(in) :: atmtab_send(:),atm_indx_send(:)
4510 : integer,intent(out),allocatable :: buf_int(:)
4511 : real(dp),intent(out),allocatable :: buf_dp(:)
4512 : type(pawrhoij_type),target,intent(in) :: pawrhoij(:)
4513 :
4514 : !Local variables-------------------------------
4515 : !scalars
4516 : integer :: cplex,ii,indx_int,indx_dp, iatom_tot,irhoij,irhoij_send,isp,jj,lmn2_size,lmnmix
4517 : integer :: ngrhoij,nselect,nspden,qphase,rhoij_size2
4518 : integer :: use_rhoijp,use_rhoijres,use_rhoij_
4519 : character(len=500) :: msg
4520 : type(pawrhoij_type),pointer :: pawrhoij1
4521 : !arrays
4522 :
4523 : ! *********************************************************************
4524 :
4525 : !Compute sizes of buffers
4526 8 : buf_int_size=0;buf_dp_size=0
4527 8 : nselect=0;lmnmix=0;ngrhoij=0;rhoij_size2=0
4528 8 : use_rhoijp=0;use_rhoijres=0;use_rhoij_=0
4529 16 : do irhoij_send=1,nrhoij_send
4530 8 : iatom_tot=atmtab_send(irhoij_send)
4531 8 : irhoij=atm_indx_send(iatom_tot)
4532 8 : if (irhoij == -1) then
4533 0 : msg="Error in pawrhoij_isendreceive_fillbuffer atom not found"
4534 0 : LIBPAW_BUG(msg)
4535 : end if
4536 8 : pawrhoij1=>pawrhoij(irhoij)
4537 8 : cplex =pawrhoij1%cplex_rhoij
4538 8 : qphase =pawrhoij1%qphase
4539 8 : lmn2_size=pawrhoij1%lmn2_size
4540 8 : nspden =pawrhoij1%nspden
4541 8 : lmnmix=pawrhoij1%lmnmix_sz
4542 8 : ngrhoij=pawrhoij1%ngrhoij
4543 8 : use_rhoijp=pawrhoij1%use_rhoijp
4544 8 : use_rhoijres=pawrhoij1%use_rhoijres
4545 8 : use_rhoij_=pawrhoij1%use_rhoij_
4546 8 : buf_int_size=buf_int_size+16
4547 8 : if (use_rhoijp>0) then
4548 8 : nselect=pawrhoij1%nrhoijsel
4549 8 : buf_int_size=buf_int_size+nselect
4550 8 : buf_dp_size=buf_dp_size + cplex*qphase*nselect*nspden
4551 : end if
4552 8 : if (lmnmix>0) buf_int_size=buf_int_size+lmnmix
4553 8 : if (ngrhoij>0) buf_dp_size=buf_dp_size + cplex*qphase*lmn2_size*nspden*ngrhoij
4554 8 : if (use_rhoijres>0) buf_dp_size=buf_dp_size + cplex*qphase*lmn2_size*nspden
4555 16 : if (use_rhoij_>0) then
4556 0 : rhoij_size2=size(pawrhoij1%rhoij_,dim=2)
4557 0 : buf_dp_size=buf_dp_size + cplex*qphase*lmn2_size*rhoij_size2
4558 : end if
4559 : end do
4560 :
4561 : !Fill input buffers
4562 24 : LIBPAW_ALLOCATE(buf_int,(buf_int_size))
4563 24 : LIBPAW_ALLOCATE(buf_dp,(buf_dp_size))
4564 8 : indx_int=1;indx_dp =1
4565 8 : lmnmix=0;ngrhoij=0;nselect=0;rhoij_size2=0
4566 8 : use_rhoijp=0;use_rhoijres=0;use_rhoij_=0
4567 16 : do irhoij_send=1,nrhoij_send
4568 8 : iatom_tot=atmtab_send(irhoij_send)
4569 8 : irhoij=atm_indx_send(iatom_tot)
4570 8 : pawrhoij1=>pawrhoij(irhoij)
4571 8 : cplex =pawrhoij1%cplex_rhoij
4572 8 : qphase =pawrhoij1%qphase
4573 8 : lmn2_size=pawrhoij1%lmn2_size
4574 8 : nspden =pawrhoij1%nspden
4575 8 : lmnmix=pawrhoij1%lmnmix_sz
4576 8 : ngrhoij=pawrhoij1%ngrhoij
4577 8 : use_rhoijp=pawrhoij1%use_rhoijp
4578 8 : nselect=pawrhoij1%nrhoijsel
4579 8 : use_rhoijres=pawrhoij1%use_rhoijres
4580 8 : use_rhoij_ =pawrhoij1%use_rhoij_
4581 8 : if (use_rhoij_ > 0) then
4582 0 : rhoij_size2 =size(pawrhoij1%rhoij_,dim=2)
4583 : end if
4584 8 : buf_int(indx_int)=atmtab_send(irhoij_send) ;indx_int=indx_int+1
4585 8 : buf_int(indx_int)=cplex ;indx_int=indx_int+1
4586 8 : buf_int(indx_int)=qphase ;indx_int=indx_int+1
4587 8 : buf_int(indx_int)=lmn2_size ;indx_int=indx_int+1
4588 8 : buf_int(indx_int)=nspden ;indx_int=indx_int+1
4589 8 : buf_int(indx_int)=nselect ;indx_int=indx_int+1
4590 8 : buf_int(indx_int)=lmnmix ;indx_int=indx_int+1
4591 8 : buf_int(indx_int)=ngrhoij ;indx_int=indx_int+1
4592 8 : buf_int(indx_int)=use_rhoijp ;indx_int=indx_int+1
4593 8 : buf_int(indx_int)=use_rhoijres ;indx_int=indx_int+1
4594 8 : buf_int(indx_int)=use_rhoij_ ;indx_int=indx_int+1
4595 8 : buf_int(indx_int)=rhoij_size2 ;indx_int=indx_int+1
4596 8 : buf_int(indx_int)=pawrhoij1%itypat ;indx_int=indx_int+1
4597 8 : buf_int(indx_int)=pawrhoij1%lmn_size ;indx_int=indx_int+1
4598 8 : buf_int(indx_int)=pawrhoij1%nsppol ;indx_int=indx_int+1
4599 8 : buf_int(indx_int)=pawrhoij1%nspinor ;indx_int=indx_int+1
4600 8 : if (use_rhoijp>0) then
4601 104 : buf_int(indx_int:indx_int+nselect-1)=pawrhoij1%rhoijselect(1:nselect)
4602 16 : indx_int=indx_int+nselect
4603 16 : do isp=1,nspden
4604 24 : do ii=1,qphase
4605 8 : jj=(ii-1)*cplex*lmn2_size
4606 104 : buf_dp(indx_dp:indx_dp+cplex*nselect-1)=pawrhoij1%rhoijp(jj+1:jj+cplex*nselect,isp)
4607 16 : indx_dp=indx_dp+cplex*nselect
4608 : end do
4609 : end do
4610 : end if
4611 8 : if (lmnmix>0) then
4612 0 : buf_int(indx_int:indx_int+lmnmix-1)=pawrhoij1%kpawmix(1:lmnmix)
4613 : indx_int=indx_int+lmnmix
4614 : end if
4615 8 : if (ngrhoij>0) then
4616 0 : do isp=1,nspden
4617 0 : do ii=1,cplex*qphase*lmn2_size
4618 0 : buf_dp(indx_dp:indx_dp+ngrhoij-1)=pawrhoij1%grhoij(1:ngrhoij,ii,isp)
4619 0 : indx_dp=indx_dp+ngrhoij
4620 : end do
4621 : end do
4622 : end if
4623 8 : if (use_rhoijres>0) then
4624 0 : do isp=1,nspden
4625 0 : buf_dp(indx_dp:indx_dp+cplex*qphase*lmn2_size-1)=pawrhoij1%rhoijres(1:cplex*qphase*lmn2_size,isp)
4626 0 : indx_dp=indx_dp+cplex*qphase*lmn2_size
4627 : end do
4628 : end if
4629 16 : if (use_rhoij_>0) then
4630 0 : do isp=1,rhoij_size2
4631 0 : buf_dp(indx_dp:indx_dp+cplex*qphase*lmn2_size-1)=pawrhoij1%rhoij_(1:cplex*qphase*lmn2_size,isp)
4632 0 : indx_dp=indx_dp+cplex*qphase*lmn2_size
4633 : end do
4634 : end if
4635 : end do !irhoij_send
4636 :
4637 : !Check
4638 8 : if ((indx_int-1/=buf_int_size).or.(indx_dp-1/=buf_dp_size)) then
4639 0 : write(msg,'(a,i10,a,i10)') 'Wrong buffer sizes: buf_int_size=',buf_int_size,' buf_dp_size=',buf_dp_size
4640 0 : LIBPAW_BUG(msg)
4641 : end if
4642 :
4643 8 : end subroutine pawrhoij_isendreceive_fillbuffer
4644 : !!***
4645 :
4646 : !----------------------------------------------------------------------
4647 :
4648 0 : END MODULE m_pawrhoij
4649 : !!***
|
