TABLE OF CONTENTS
- ABINIT/m_pawrhoij
- m_pawdij/pawrhoij_print_rhoij
- m_pawrhoij/pawrhoij_alloc
- m_pawrhoij/pawrhoij_bcast
- m_pawrhoij/pawrhoij_copy
- m_pawrhoij/pawrhoij_filter
- m_pawrhoij/pawrhoij_free
- m_pawrhoij/pawrhoij_free_unpacked
- m_pawrhoij/pawrhoij_gather
- m_pawrhoij/pawrhoij_init_unpacked
- m_pawrhoij/pawrhoij_inquire_dim
- m_pawrhoij/pawrhoij_io
- m_pawrhoij/pawrhoij_isendreceive_fillbuffer
- m_pawrhoij/pawrhoij_isendreceive_getbuffer
- m_pawrhoij/pawrhoij_mpisum_unpacked_1D
- m_pawrhoij/pawrhoij_mpisum_unpacked_2D
- m_pawrhoij/pawrhoij_nullify
- m_pawrhoij/pawrhoij_redistribute
- m_pawrhoij/pawrhoij_symrhoij
- m_pawrhoij/pawrhoij_type
- m_pawrhoij/pawrhoij_unpack
ABINIT/m_pawrhoij [ Modules ]
NAME
m_pawrhoij
FUNCTION
This module contains the definition of the pawrhoij_type structured datatype, as well as related functions and methods. pawrhoij_type variables define rhoij occupancies matrixes used within PAW formalism.
COPYRIGHT
Copyright (C) 2012-2024 ABINIT group (MT, FJ) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt .
NOTES
FOR DEVELOPPERS: in order to preserve the portability of libPAW library, please consult ~abinit/src/??_libpaw/libpaw-coding-rules.txt
SOURCE
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
m_pawdij/pawrhoij_print_rhoij [ Functions ]
[ Top ] [ m_pawdij ] [ Functions ]
NAME
pawrhoij_print_rhoij
FUNCTION
Print out the content of a Rho_ij matrix (occupation matrix) in a suitable format
INPUTS
rhoij(cplex*lmn2_size,nspden)= input matrix to be printed cplex=1 if Rhoij is real, 2 if Rhoij is complex qphase=2 if rhoij has a exp(iqR) phase, 1 if not iatom=current atom natom=total number of atoms in the system [opt_prtvol]= >=0 if up to 12 components of _ij matrix have to be printed <0 if all components of ij_ matrix have to be printed (optional) [mode_paral]= parallel printing mode (optional, default='COLL') [test_value]=(real number) if positive, print a warning when the magnitude of Dij is greater (optional) [l_only]=if >=0 only parts of rhoij corresponding to li=l_only are printed (optional); Needs indlmn(:,:) optional argument. [title_msg]=message to print as title (optional) [unit]=the unit number for output (optional) [indlmn(6,lmn_size)]= array giving l,m,n,lm,ln,s
OUTPUT
(Only writing)
NOTES
SOURCE
3285 subroutine pawrhoij_print_rhoij(rhoij,cplex,qphase,iatom,natom,& 3286 & rhoijselect,test_value,title_msg,unit,opt_prtvol,& 3287 & 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 integer,pointer :: l_index(:),my_rhoijselect(:) 3312 real(dp),pointer :: rhoij_(:) 3313 3314 ! ************************************************************************* 3315 3316 !Optional arguments 3317 my_unt =std_out ; if (PRESENT(unit )) my_unt =unit 3318 my_mode ='COLL' ; if (PRESENT(mode_paral)) my_mode =mode_paral 3319 my_prtvol=1 ; if (PRESENT(opt_prtvol)) my_prtvol=opt_prtvol 3320 my_test_value=-one; if (PRESENT(test_value)) my_test_value=test_value 3321 my_l_only=-1 ; if (PRESENT(l_only)) my_l_only=l_only 3322 3323 if (my_l_only>=0.and.(.not.present(indlmn))) then 3324 msg='pawrhoij_print_rhoij: l_only>=0 and indlmn not present!' 3325 LIBPAW_BUG(msg) 3326 end if 3327 3328 !Title 3329 if (present(title_msg)) then 3330 if (trim(title_msg)/='') then 3331 write(msg, '(2a)') ch10,trim(title_msg) 3332 call wrtout(my_unt,msg,my_mode) 3333 end if 3334 end if 3335 3336 !Inits 3337 my_nspden=size(rhoij,2) 3338 my_lmn2_size=size(rhoij,1)/cplex/qphase 3339 my_lmn_size=int(dsqrt(two*dble(my_lmn2_size))) 3340 my_cplex=merge(cplex,2,qphase==1) 3341 my_opt_sym=2 3342 3343 !Packed storage 3344 my_opt_pack=0 3345 rhoij_size=my_lmn2_size 3346 my_rhoijselect => idum 3347 if (PRESENT(rhoijselect)) then 3348 nrhoijsel=count(rhoijselect(:)>0) 3349 if (nrhoijsel>0) then 3350 my_opt_pack=1 3351 rhoij_size=nrhoijsel 3352 my_rhoijselect => rhoijselect(1:nrhoijsel) 3353 end if 3354 end if 3355 3356 if (my_l_only<0) then 3357 l_index => idum 3358 else 3359 LIBPAW_POINTER_ALLOCATE(l_index,(my_lmn_size)) 3360 do kk=1,my_lmn_size 3361 l_index(kk)=indlmn(1,kk) 3362 end do 3363 end if 3364 3365 if (qphase==2) then 3366 LIBPAW_DATATYPE_ALLOCATE(rhoij_,(2*rhoij_size)) 3367 end if 3368 3369 ! === Loop over Rho_ij components === 3370 do irhoij=1,my_nspden 3371 3372 !Rebuild rhoij according to qphase 3373 if (qphase==1) then 3374 rhoij_ => rhoij(1:cplex*rhoij_size,irhoij) 3375 else 3376 if (cplex==1) then 3377 do kk=1,rhoij_size 3378 rhoij_(2*kk-1)=rhoij(kk,irhoij) 3379 rhoij_(2*kk )=rhoij(kk+my_lmn2_size,irhoij) 3380 end do 3381 else 3382 do kk=1,rhoij_size 3383 rhoij_(2*kk-1)=rhoij(2*kk-1,irhoij)-rhoij(2*kk +2*my_lmn2_size,irhoij) 3384 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 if (natom>1.or.my_nspden>1) then 3391 if (my_l_only<0) then 3392 if (my_nspden==1) write(msg,'(a,i3)') ' Atom #',iatom 3393 if (my_nspden==2) write(msg,'(a,i3,a,i1)')' Atom #',iatom,' - Spin component ',irhoij 3394 if (my_nspden==4) write(msg,'(a,i3,2a)') ' Atom #',iatom,' - Component ',trim(dspin(irhoij+2*(my_nspden/4))) 3395 else 3396 if (my_nspden==1) write(msg,'(a,i3,a,i1,a)') ' Atom #',iatom,& 3397 & ' - L=',my_l_only,' ONLY' 3398 if (my_nspden==2) write(msg,'(a,i3,a,i1,a,i1)')' Atom #',iatom,& 3399 & ' - L=',my_l_only,' ONLY - Spin component ',irhoij 3400 if (my_nspden==4) write(msg,'(a,i3,a,i1,3a)') ' Atom #',iatom,& 3401 & ' - L=',my_l_only,' ONLY - Component ',trim(dspin(irhoij+2*(my_nspden/4))) 3402 end if 3403 call wrtout(my_unt,msg,my_mode) 3404 else if (my_l_only>=0) then 3405 write(msg,'(a,i1,a)') ' L=',my_l_only,' ONLY' 3406 call wrtout(my_unt,msg,my_mode) 3407 end if 3408 3409 !Printing 3410 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 & mode_paral=my_mode) 3414 3415 end do !irhoij 3416 3417 if (qphase==2) then 3418 LIBPAW_DATATYPE_DEALLOCATE(rhoij_) 3419 end if 3420 if (my_l_only>=0) then 3421 LIBPAW_POINTER_DEALLOCATE(l_index) 3422 end if 3423 3424 end subroutine pawrhoij_print_rhoij
m_pawrhoij/pawrhoij_alloc [ Functions ]
[ Top ] [ m_pawrhoij ] [ Functions ]
NAME
pawrhoij_alloc
FUNCTION
Initialize and allocate a pawrhoij datastructure
INPUTS
[comm_atom] = communicator over atoms (OPTIONAL)
cplex_rhoij=1 if rhoij are real, 2 if rhoij are complex (spin-orbit, pawcpxocc=2, ...)
[my_atmtab(:)] = Index of atoms treated by current proc (OPTIONAL)
nspden=number of spin-components for rhoij
nsppol=number of spinorial components for rhoij
nsppol=number of independant spin-components for rhoij
typat(:)=types of atoms
[lmnsize(:)]=array of (l,m,n) sizes for rhoij for each type of atom (OPTIONAL)
must be present if [pawtab] argument is not passed
[ngrhoij]=number of gradients to be allocated (OPTIONAL, default=0)
[nlmnmix]=number of rhoij elements to be mixed during SCF cycle (OPTIONAL, default=0)
[pawtab(:)] <type(pawtab_type)>=paw tabulated starting data (OPTIONAL)
must be present if [lmnsize(:)] argument is not passed
[qphase]=2 if the rhoij contain a exp(iqR) phase, 1 otherwise (OPTIONAL, default=1)
Typical use: 1st-order rhoij at q<>0
[use_rhoij_]=1 if pawrhoij(:)%rhoij_ has to be allocated (OPTIONAL, default=0)
[use_rhoijp]=1 if pawrhoij(:)%rhoijp has to be allocated (OPTIONAL, default=1)
(in that case, pawrhoij%rhoijselect is also allocated)
[use_rhoijres]=1 if pawrhoij(:)%rhoijres has to be allocated (OPTIONAL, default=0)
SIDE EFFECTS
pawrhoij(:)<type(pawrhoij_type)>= rhoij datastructure
NOTES
One of the two optional arguments lmnsize(:) or pawtab(:) must be present ! If both are present, only pawtab(:) is used.
SOURCE
235 subroutine pawrhoij_alloc(pawrhoij,cplex_rhoij,nspden,nspinor,nsppol,typat,& 236 & lmnsize,ngrhoij,nlmnmix,pawtab,qphase,use_rhoij_,use_rhoijp,& ! Optional 237 & 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 integer,pointer :: lmn_size(:),my_atmtab(:) 258 259 ! ************************************************************************* 260 261 nrhoij=size(pawrhoij);natom=size(typat) 262 if (nrhoij>natom) then 263 msg=' wrong sizes (1) !' 264 LIBPAW_BUG(msg) 265 end if 266 267 !Select lmn_size for each atom type 268 if (present(pawtab)) then 269 nn1=size(pawtab) 270 if (maxval(typat)>nn1) then 271 msg=' wrong sizes (2) !' 272 LIBPAW_BUG(msg) 273 end if 274 LIBPAW_POINTER_ALLOCATE(lmn_size,(nn1)) 275 do itypat=1,nn1 276 lmn_size(itypat)=pawtab(itypat)%lmn_size 277 end do 278 else if (present(lmnsize)) then 279 nn1=size(lmnsize) 280 if (maxval(typat)>nn1) then 281 msg=' wrong sizes (3) !' 282 LIBPAW_BUG(msg) 283 end if 284 lmn_size => lmnsize 285 else 286 msg=' one of the 2 arguments pawtab or lmnsize must be present !' 287 LIBPAW_BUG(msg) 288 end if 289 290 !Set up parallelism over atoms 291 paral_atom=(present(comm_atom).and.(nrhoij/=natom)) 292 nullify(my_atmtab);if (present(mpi_atmtab)) my_atmtab => mpi_atmtab 293 my_comm_atom=xmpi_comm_self;if (present(comm_atom)) my_comm_atom=comm_atom 294 call get_my_atmtab(my_comm_atom,my_atmtab,my_atmtab_allocated,paral_atom,natom) 295 296 my_qphase=1;if (present(qphase)) my_qphase=qphase 297 298 if (nrhoij>0) then 299 do irhoij=1,nrhoij 300 irhoij_=irhoij;if (paral_atom) irhoij_=my_atmtab(irhoij) 301 itypat=typat(irhoij_) 302 303 lmn2_size=lmn_size(itypat)*(lmn_size(itypat)+1)/2 304 305 ! Scalars initializations 306 pawrhoij(irhoij)%cplex_rhoij=cplex_rhoij 307 pawrhoij(irhoij)%qphase=my_qphase 308 pawrhoij(irhoij)%itypat=itypat 309 pawrhoij(irhoij)%lmn_size=lmn_size(itypat) 310 pawrhoij(irhoij)%lmn2_size=lmn2_size 311 pawrhoij(irhoij)%nspden=nspden 312 pawrhoij(irhoij)%nspinor=nspinor 313 pawrhoij(irhoij)%nsppol=nsppol 314 pawrhoij(irhoij)%nrhoijsel=0 315 pawrhoij(irhoij)%lmnmix_sz=0 316 pawrhoij(irhoij)%ngrhoij=0 317 pawrhoij(irhoij)%use_rhoij_=0 318 pawrhoij(irhoij)%use_rhoijres=0 319 320 ! Arrays allocations 321 has_rhoijp=.true.; if (present(use_rhoijp)) has_rhoijp=(use_rhoijp>0) 322 if (has_rhoijp) then 323 pawrhoij(irhoij)%use_rhoijp=1 324 LIBPAW_ALLOCATE(pawrhoij(irhoij)%rhoijselect,(lmn2_size)) 325 LIBPAW_ALLOCATE(pawrhoij(irhoij)%rhoijp,(cplex_rhoij*my_qphase*lmn2_size,nspden)) 326 pawrhoij(irhoij)%rhoijselect(:)=0 327 pawrhoij(irhoij)%rhoijp(:,:)=zero 328 end if 329 330 if (present(ngrhoij)) then 331 if (ngrhoij>0) then 332 pawrhoij(irhoij)%ngrhoij=ngrhoij 333 LIBPAW_ALLOCATE(pawrhoij(irhoij)%grhoij,(ngrhoij,cplex_rhoij*my_qphase*lmn2_size,nspden)) 334 pawrhoij(irhoij)%grhoij=zero 335 end if 336 end if 337 if (present(nlmnmix)) then 338 if (nlmnmix>0) then 339 pawrhoij(irhoij)%lmnmix_sz=nlmnmix 340 LIBPAW_ALLOCATE(pawrhoij(irhoij)%kpawmix,(nlmnmix)) 341 pawrhoij(irhoij)%kpawmix=0 342 end if 343 end if 344 if (present(use_rhoij_)) then 345 if (use_rhoij_>0) then 346 pawrhoij(irhoij)%use_rhoij_=use_rhoij_ 347 LIBPAW_ALLOCATE(pawrhoij(irhoij)%rhoij_,(cplex_rhoij*my_qphase*lmn2_size,nspden)) 348 pawrhoij(irhoij)%rhoij_=zero 349 end if 350 end if 351 if (present(use_rhoijres)) then 352 if (use_rhoijres>0) then 353 pawrhoij(irhoij)%use_rhoijres=use_rhoijres 354 LIBPAW_ALLOCATE(pawrhoij(irhoij)%rhoijres,(cplex_rhoij*my_qphase*lmn2_size,nspden)) 355 pawrhoij(irhoij)%rhoijres=zero 356 end if 357 end if 358 359 end do 360 end if 361 362 if (present(pawtab)) then 363 LIBPAW_POINTER_DEALLOCATE(lmn_size) 364 end if 365 366 !Destroy atom table used for parallelism 367 call free_my_atmtab(my_atmtab,my_atmtab_allocated) 368 369 end subroutine pawrhoij_alloc
m_pawrhoij/pawrhoij_bcast [ Functions ]
[ Top ] [ m_pawrhoij ] [ Functions ]
NAME
pawrhoij_bcast
FUNCTION
Broadcast pawrhoij datastructures Can take into account a distribution of data over a "atom" communicator
INPUTS
master=master communicator receiving data mpicomm= MPI communicator comm_atom= --optional-- MPI communicator over atoms pawrhoij_in(:)<type(pawrhoij_type)>= input rhoij datastructures on master process
OUTPUT
pawrhoij_out(:)<type(pawrhoij_type)>= output rhoij datastructure on every process Eventually distributed according to comm_atom communicator
SOURCE
1541 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 integer,allocatable :: atmtab(:),buf_int_size_i(:),buf_dp_size_i(:) 1561 integer,allocatable :: count_dp(:),count_int(:),disp_dp(:),disp_int(:),nrhoij_out_i(:) 1562 integer,allocatable,target :: buf_int(:) 1563 integer,pointer :: buf_int_all(:),my_atmtab(:) 1564 real(dp),allocatable,target :: buf_dp(:) 1565 real(dp),pointer :: buf_dp_all(:) 1566 1567 ! ************************************************************************* 1568 1569 !Load MPI "atom" distribution data 1570 my_comm_atom=xmpi_comm_self;nproc_atom=1;me_atom=0 1571 if (present(comm_atom)) then 1572 my_comm_atom=comm_atom 1573 me_atom=xmpi_comm_rank(my_comm_atom) 1574 nproc_atom=xmpi_comm_size(my_comm_atom) 1575 paral_atom=(nproc_atom>1) 1576 if (my_comm_atom/=mpicomm.and.nproc_atom/=1) then 1577 msg='wrong comm_atom communicator !' 1578 LIBPAW_BUG(msg) 1579 end if 1580 end if 1581 1582 !Load global MPI data 1583 me=xmpi_comm_rank(mpicomm) 1584 nproc=xmpi_comm_size(mpicomm) 1585 1586 !Just copy in case of a sequential run 1587 if (nproc==1.and.nproc_atom==1) then 1588 call pawrhoij_copy(pawrhoij_in,pawrhoij_out,.false.,.false.,.false.) 1589 return 1590 end if 1591 1592 !Retrieve and test pawrhoij sizes 1593 nrhoij_in=0;if (me==master) nrhoij_in=size(pawrhoij_in) 1594 nrhoij_out=size(pawrhoij_out);nrhoij_out_all=nrhoij_out 1595 if (paral_atom) then 1596 LIBPAW_ALLOCATE(nrhoij_out_i,(nproc_atom)) 1597 buf_size(1)=nrhoij_out 1598 call xmpi_allgather(buf_size,1,nrhoij_out_i,my_comm_atom,ierr) 1599 nrhoij_out_all=sum(nrhoij_out_i) 1600 end if 1601 if (me==master.and.nrhoij_in/=nrhoij_out_all) then 1602 msg='pawrhoij_in or pawrhoij_out wrongly allocated!' 1603 LIBPAW_BUG(msg) 1604 end if 1605 1606 !Retrieve table(s) of atoms (if necessary) 1607 if (paral_atom) then 1608 nullify(my_atmtab) 1609 call get_my_atmtab(my_comm_atom,my_atmtab,my_atmtab_allocated,paral_atom, & 1610 & nrhoij_out_all,my_natom_ref=nrhoij_out) 1611 LIBPAW_ALLOCATE(disp_int,(nproc_atom)) 1612 disp_int(1)=0 1613 do iproc=2,nproc_atom 1614 disp_int(iproc)=disp_int(iproc-1)+nrhoij_out_i(iproc-1) 1615 end do 1616 LIBPAW_ALLOCATE(atmtab,(nrhoij_in)) 1617 call xmpi_gatherv(my_atmtab,nrhoij_out,atmtab,nrhoij_out_i,disp_int,& 1618 & master,my_comm_atom,ierr) 1619 LIBPAW_DEALLOCATE(disp_int) 1620 end if 1621 1622 !Compute sizes of input buffers and broadcast them 1623 LIBPAW_ALLOCATE(buf_int_size_i,(nrhoij_out_all)) 1624 LIBPAW_ALLOCATE(buf_dp_size_i ,(nrhoij_out_all)) 1625 if (me==master) then 1626 buf_int_size_i(:)=0;buf_dp_size_i(:)=0 1627 do irhoij=1,nrhoij_in 1628 jrhoij=irhoij;if (paral_atom) jrhoij=atmtab(irhoij) 1629 cplex =pawrhoij_in(jrhoij)%cplex_rhoij 1630 qphase =pawrhoij_in(jrhoij)%qphase 1631 lmn2_size =pawrhoij_in(jrhoij)%lmn2_size 1632 nspden =pawrhoij_in(jrhoij)%nspden 1633 lmnmix =pawrhoij_in(jrhoij)%lmnmix_sz 1634 ngrhoij =pawrhoij_in(jrhoij)%ngrhoij 1635 use_rhoijp =pawrhoij_in(jrhoij)%use_rhoijp 1636 use_rhoijres=pawrhoij_in(jrhoij)%use_rhoijres 1637 use_rhoij_ =pawrhoij_in(jrhoij)%use_rhoij_ 1638 buf_int_size_i(irhoij)=buf_int_size_i(irhoij)+16 1639 if (ngrhoij>0) buf_dp_size_i(irhoij)=buf_dp_size_i(irhoij)+cplex*qphase*lmn2_size*nspden*ngrhoij 1640 if (use_rhoijres>0) buf_dp_size_i(irhoij)=buf_dp_size_i(irhoij)+cplex*qphase*lmn2_size*nspden 1641 if (use_rhoijp>0) then 1642 nselect=pawrhoij_in(jrhoij)%nrhoijsel 1643 buf_int_size_i(irhoij)=buf_int_size_i(irhoij)+nselect 1644 buf_dp_size_i(irhoij)=buf_dp_size_i(irhoij)+cplex*qphase*nselect*nspden 1645 end if 1646 if (use_rhoij_>0) then 1647 rhoij_size2=size(pawrhoij_in(jrhoij)%rhoij_,dim=2) 1648 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 call xmpi_bcast(buf_int_size_i,master,mpicomm,ierr) 1653 call xmpi_bcast(buf_dp_size_i,master,mpicomm,ierr) 1654 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 if (paral_atom) then 1658 LIBPAW_ALLOCATE(count_int,(nproc_atom)) 1659 LIBPAW_ALLOCATE(count_dp,(nproc_atom)) 1660 LIBPAW_ALLOCATE(disp_int,(nproc_atom)) 1661 LIBPAW_ALLOCATE(disp_dp,(nproc_atom)) 1662 indx_int=0 1663 do iproc=1,nproc_atom 1664 ii=nrhoij_out_i(iproc) 1665 count_int(iproc)=sum(buf_int_size_i(indx_int+1:indx_int+ii)) 1666 count_dp (iproc)=sum(buf_dp_size_i (indx_int+1:indx_int+ii)) 1667 indx_int=indx_int+ii 1668 end do 1669 disp_int(1)=0;disp_dp(1)=0 1670 do iproc=2,nproc_atom 1671 disp_int(iproc)=disp_int(iproc-1)+count_int(iproc-1) 1672 disp_dp (iproc)=disp_dp (iproc-1)+count_dp (iproc-1) 1673 end do 1674 if (buf_int_size_all/=sum(count_int).or.buf_dp_size_all/=sum(count_dp)) then 1675 msg='(1) Wrong buffer sizes !' 1676 LIBPAW_BUG(msg) 1677 end if 1678 buf_int_size=count_int(me_atom+1) 1679 buf_dp_size =count_dp(me_atom+1) 1680 LIBPAW_ALLOCATE(buf_int,(buf_int_size)) 1681 LIBPAW_ALLOCATE(buf_dp ,(buf_dp_size)) 1682 if (me==master) then 1683 LIBPAW_POINTER_ALLOCATE(buf_int_all,(buf_int_size_all)) 1684 LIBPAW_POINTER_ALLOCATE(buf_dp_all ,(buf_dp_size_all)) 1685 else 1686 LIBPAW_POINTER_ALLOCATE(buf_int_all,(1)) 1687 LIBPAW_POINTER_ALLOCATE(buf_dp_all ,(1)) 1688 end if 1689 LIBPAW_DEALLOCATE(nrhoij_out_i) 1690 else 1691 buf_int_size=buf_int_size_all 1692 buf_dp_size =buf_dp_size_all 1693 LIBPAW_ALLOCATE(buf_int,(buf_int_size)) 1694 LIBPAW_ALLOCATE(buf_dp ,(buf_dp_size)) 1695 buf_int_all => buf_int 1696 buf_dp_all => buf_dp 1697 end if 1698 LIBPAW_DEALLOCATE(buf_int_size_i) 1699 LIBPAW_DEALLOCATE(buf_dp_size_i) 1700 1701 !Fill input buffers 1702 if (me==master) then 1703 indx_int=1;indx_dp =1 1704 do irhoij=1,nrhoij_in 1705 jrhoij=irhoij;if (paral_atom) jrhoij=atmtab(irhoij) 1706 cplex =pawrhoij_in(jrhoij)%cplex_rhoij 1707 qphase =pawrhoij_in(jrhoij)%qphase 1708 lmn2_size =pawrhoij_in(jrhoij)%lmn2_size 1709 nspden =pawrhoij_in(jrhoij)%nspden 1710 lmnmix =pawrhoij_in(jrhoij)%lmnmix_sz 1711 ngrhoij =pawrhoij_in(jrhoij)%ngrhoij 1712 use_rhoijp =pawrhoij_in(jrhoij)%use_rhoijp 1713 nselect =pawrhoij_in(jrhoij)%nrhoijsel 1714 use_rhoijres=pawrhoij_in(jrhoij)%use_rhoijres 1715 use_rhoij_ =pawrhoij_in(jrhoij)%use_rhoij_ 1716 rhoij_size2 =size(pawrhoij_in(jrhoij)%rhoij_,dim=2) 1717 buf_int_all(indx_int)=jrhoij ;indx_int=indx_int+1 ! Not used ! 1718 buf_int_all(indx_int)=cplex ;indx_int=indx_int+1 1719 buf_int_all(indx_int)=qphase ;indx_int=indx_int+1 1720 buf_int_all(indx_int)=lmn2_size ;indx_int=indx_int+1 1721 buf_int_all(indx_int)=nspden ;indx_int=indx_int+1 1722 buf_int_all(indx_int)=nselect ;indx_int=indx_int+1 1723 buf_int_all(indx_int)=lmnmix ;indx_int=indx_int+1 1724 buf_int_all(indx_int)=ngrhoij ;indx_int=indx_int+1 1725 buf_int_all(indx_int)=use_rhoijp ;indx_int=indx_int+1 1726 buf_int_all(indx_int)=use_rhoijres ;indx_int=indx_int+1 1727 buf_int_all(indx_int)=use_rhoij_ ;indx_int=indx_int+1 1728 buf_int_all(indx_int)=rhoij_size2 ;indx_int=indx_int+1 1729 buf_int_all(indx_int)=pawrhoij_in(jrhoij)%itypat ;indx_int=indx_int+1 1730 buf_int_all(indx_int)=pawrhoij_in(jrhoij)%lmn_size;indx_int=indx_int+1 1731 buf_int_all(indx_int)=pawrhoij_in(jrhoij)%nsppol ;indx_int=indx_int+1 1732 buf_int_all(indx_int)=pawrhoij_in(jrhoij)%nspinor ;indx_int=indx_int+1 1733 if (use_rhoijp>0) then 1734 buf_int_all(indx_int:indx_int+nselect-1)=pawrhoij_in(jrhoij)%rhoijselect(1:nselect) 1735 indx_int=indx_int+nselect 1736 do isp=1,nspden 1737 do ii=1,qphase 1738 jj=(ii-1)*cplex*lmn2_size 1739 buf_dp_all(indx_dp:indx_dp+cplex*nselect-1)= & 1740 & pawrhoij_in(jrhoij)%rhoijp(jj+1:jj+cplex*nselect,isp) 1741 indx_dp=indx_dp+cplex*nselect 1742 end do 1743 end do 1744 end if 1745 if (lmnmix>0) then 1746 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 if (ngrhoij>0) then 1750 do isp=1,nspden 1751 do ii=1,cplex*qphase*lmn2_size 1752 buf_dp_all(indx_dp:indx_dp+ngrhoij-1)=pawrhoij_in(jrhoij)%grhoij(1:ngrhoij,ii,isp) 1753 indx_dp=indx_dp+ngrhoij 1754 end do 1755 end do 1756 end if 1757 if (use_rhoijres>0) then 1758 do isp=1,nspden 1759 buf_dp_all(indx_dp:indx_dp+cplex*qphase*lmn2_size-1)= & 1760 & pawrhoij_in(jrhoij)%rhoijres(1:cplex*qphase*lmn2_size,isp) 1761 indx_dp=indx_dp+cplex*qphase*lmn2_size 1762 end do 1763 end if 1764 if (use_rhoij_>0) then 1765 do isp=1,rhoij_size2 1766 buf_dp_all(indx_dp:indx_dp+cplex*qphase*lmn2_size-1)= & 1767 & pawrhoij_in(jrhoij)%rhoij_(1:cplex*qphase*lmn2_size,isp) 1768 indx_dp=indx_dp+cplex*qphase*lmn2_size 1769 end do 1770 end if 1771 end do 1772 ! Check 1773 if ((indx_int-1/=buf_int_size_all).or.(indx_dp-1/=buf_dp_size_all)) then 1774 msg='(2) Wrong buffer sizes !' 1775 LIBPAW_BUG(msg) 1776 end if 1777 end if ! me=master 1778 1779 !Communicate 1780 if (paral_atom) then 1781 call xmpi_scatterv(buf_int_all,count_int,disp_int,buf_int,buf_int_size,master,mpicomm,ierr) 1782 call xmpi_scatterv(buf_dp_all ,count_dp ,disp_dp ,buf_dp ,buf_dp_size ,master,mpicomm,ierr) 1783 else 1784 call xmpi_bcast(buf_int,master,mpicomm,ierr) 1785 call xmpi_bcast(buf_dp ,master,mpicomm,ierr) 1786 end if 1787 1788 !Retrieve data from output buffer 1789 indx_int=1;indx_dp=1 1790 call pawrhoij_free(pawrhoij_out) 1791 do irhoij=1,nrhoij_out 1792 jrhoij =buf_int(indx_int);indx_int=indx_int+1 ! Not used ! 1793 cplex =buf_int(indx_int);indx_int=indx_int+1 1794 qphase =buf_int(indx_int);indx_int=indx_int+1 1795 lmn2_size =buf_int(indx_int);indx_int=indx_int+1 1796 nspden =buf_int(indx_int);indx_int=indx_int+1 1797 nselect =buf_int(indx_int);indx_int=indx_int+1 1798 lmnmix =buf_int(indx_int);indx_int=indx_int+1 1799 ngrhoij =buf_int(indx_int);indx_int=indx_int+1 1800 use_rhoijp =buf_int(indx_int);indx_int=indx_int+1 1801 use_rhoijres=buf_int(indx_int);indx_int=indx_int+1 1802 use_rhoij_ =buf_int(indx_int);indx_int=indx_int+1 1803 rhoij_size2 =buf_int(indx_int);indx_int=indx_int+1 1804 pawrhoij_out(irhoij)%itypat=buf_int(indx_int) ;indx_int=indx_int+1 1805 pawrhoij_out(irhoij)%lmn_size=buf_int(indx_int);indx_int=indx_int+1 1806 pawrhoij_out(irhoij)%nsppol=buf_int(indx_int) ;indx_int=indx_int+1 1807 pawrhoij_out(irhoij)%nspinor=buf_int(indx_int) ;indx_int=indx_int+1 1808 pawrhoij_out(irhoij)%cplex_rhoij=cplex 1809 pawrhoij_out(irhoij)%qphase=qphase 1810 pawrhoij_out(irhoij)%lmn2_size=lmn2_size 1811 pawrhoij_out(irhoij)%nspden=nspden 1812 pawrhoij_out(irhoij)%nrhoijsel=nselect 1813 pawrhoij_out(irhoij)%lmnmix_sz=lmnmix 1814 pawrhoij_out(irhoij)%ngrhoij=ngrhoij 1815 pawrhoij_out(irhoij)%use_rhoijp=use_rhoijp 1816 pawrhoij_out(irhoij)%use_rhoijres=use_rhoijres 1817 pawrhoij_out(irhoij)%use_rhoij_=use_rhoij_ 1818 if (use_rhoijp>0) then 1819 LIBPAW_ALLOCATE(pawrhoij_out(irhoij)%rhoijselect,(nselect)) 1820 pawrhoij_out(irhoij)%rhoijselect=0 1821 pawrhoij_out(irhoij)%rhoijselect(1:nselect)=buf_int(indx_int:indx_int+nselect-1) 1822 indx_int=indx_int+nselect 1823 LIBPAW_ALLOCATE(pawrhoij_out(irhoij)%rhoijp,(qphase*cplex*nselect,nspden)) 1824 do isp=1,nspden 1825 do ii=1,qphase 1826 jj=(ii-1)*cplex*lmn2_size 1827 pawrhoij_out(irhoij)%rhoijp(jj+1:jj+cplex*nselect,isp)= & 1828 & buf_dp(indx_dp:indx_dp+cplex*nselect-1) 1829 indx_dp=indx_dp+cplex*nselect 1830 end do 1831 end do 1832 end if 1833 if (lmnmix>0) then 1834 LIBPAW_ALLOCATE(pawrhoij_out(irhoij)%kpawmix,(lmnmix)) 1835 pawrhoij_out(irhoij)%kpawmix(1:lmnmix)=buf_int(indx_int:indx_int+lmnmix-1) 1836 indx_int=indx_int+lmnmix 1837 end if 1838 if (ngrhoij>0) then 1839 LIBPAW_ALLOCATE(pawrhoij_out(irhoij)%grhoij,(ngrhoij,cplex*qphase*lmn2_size,nspden)) 1840 do isp=1,nspden 1841 do ii=1,cplex*qphase*lmn2_size 1842 pawrhoij_out(irhoij)%grhoij(1:ngrhoij,ii,isp)=buf_dp(indx_dp:indx_dp+ngrhoij-1) 1843 indx_dp=indx_dp+ngrhoij 1844 end do 1845 end do 1846 end if 1847 if (use_rhoijres>0) then 1848 LIBPAW_ALLOCATE(pawrhoij_out(irhoij)%rhoijres,(cplex*qphase*lmn2_size,nspden)) 1849 do isp=1,nspden 1850 pawrhoij_out(irhoij)%rhoijres(1:cplex*qphase*lmn2_size,isp)= & 1851 & buf_dp(indx_dp:indx_dp+cplex*qphase*lmn2_size-1) 1852 indx_dp=indx_dp+cplex*qphase*lmn2_size 1853 end do 1854 end if 1855 if (use_rhoij_>0) then 1856 LIBPAW_ALLOCATE(pawrhoij_out(irhoij)%rhoij_,(cplex*qphase*lmn2_size,rhoij_size2)) 1857 do isp=1,rhoij_size2 1858 pawrhoij_out(irhoij)%rhoij_(1:cplex*qphase*lmn2_size,isp)= & 1859 & buf_dp(indx_dp:indx_dp+cplex*qphase*lmn2_size-1) 1860 indx_dp=indx_dp+cplex*qphase*lmn2_size 1861 end do 1862 end if 1863 end do 1864 !Check 1865 if ((indx_int/=1+buf_int_size).or.(indx_dp/=1+buf_dp_size)) then 1866 msg='(3) Wrong buffer sizes !' 1867 LIBPAW_BUG(msg) 1868 end if 1869 1870 !Free memory 1871 LIBPAW_DEALLOCATE(buf_int) 1872 LIBPAW_DEALLOCATE(buf_dp) 1873 if (paral_atom) then 1874 LIBPAW_POINTER_DEALLOCATE(buf_int_all) 1875 LIBPAW_POINTER_DEALLOCATE(buf_dp_all) 1876 LIBPAW_DEALLOCATE(count_int) 1877 LIBPAW_DEALLOCATE(count_dp) 1878 LIBPAW_DEALLOCATE(disp_int) 1879 LIBPAW_DEALLOCATE(disp_dp) 1880 LIBPAW_DEALLOCATE(atmtab) 1881 call free_my_atmtab(my_atmtab,my_atmtab_allocated) 1882 end if 1883 1884 end subroutine pawrhoij_bcast
m_pawrhoij/pawrhoij_copy [ Functions ]
[ Top ] [ m_pawrhoij ] [ Functions ]
NAME
pawrhoij_copy
FUNCTION
Copy one pawrhoij datastructure into another Can take into accound changes of dimensions Can copy a shared pawrhoij into distributed ones (when parallelism is activated)
INPUTS
keep_cplex= optional argument (logical, default=.TRUE.)
if .TRUE. pawrhoij_out(:)%cplex_rhoij is NOT MODIFIED,
even if different from pawrhoij_in(:)%cplex_rhoij
keep_qphase= optional argument (logical, default=.TRUE.)
if .TRUE. pawrhoij_out(:)%cplex_rhoij is NOT MODIFIED,
even if different from pawrhoij_in(:)%qphase
keep_itypat= optional argument (logical, default=.FALSE.)
if .TRUE. pawrhoij_out(:)%ityp is NOT MODIFIED,
even if different from pawrhoij_in(:)%ityp
keep_nspden= optional argument (logical, default=.TRUE.)
if .TRUE. pawrhoij_out(:)%nspden is NOT MODIFIED,
even if different from pawrhoij_in(:)%nspden
mpi_atmtab(:)=--optional-- indexes of the atoms treated by current proc
comm_atom=--optional-- MPI communicator over atoms
pawrhoij_in(:)<type(pawrhoij_type)>= input rhoij datastructure
SIDE EFFECTS
pawrhoij_out(:)<type(pawrhoij_type)>= output rhoij datastructure
NOTES
In case of a single copy operation pawrhoij_out must have been allocated.
SOURCE
518 subroutine pawrhoij_copy(pawrhoij_in,pawrhoij_cpy, & 519 & keep_cplex,keep_qphase,keep_itypat,keep_nspden,& ! optional arguments 520 & 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 integer,pointer :: my_atmtab(:) 542 integer,allocatable :: nlmn(:),typat(:) 543 type(pawrhoij_type),pointer :: pawrhoij_out(:) 544 545 ! ************************************************************************* 546 547 !Retrieve sizes 548 nrhoij_in=size(pawrhoij_in);nrhoij_out=size(pawrhoij_cpy) 549 550 !Init flags 551 keep_cplex_=.true. 552 if (present(keep_cplex)) keep_cplex_=keep_cplex 553 keep_qphase_=.true. 554 if (present(keep_qphase)) keep_qphase_=keep_qphase 555 keep_itypat_=.false. 556 if (present(keep_itypat)) keep_itypat_=keep_itypat 557 keep_nspden_=.true. 558 if (present(keep_nspden)) keep_nspden_=keep_nspden 559 560 !Set up parallelism over atoms 561 paral_atom=(present(comm_atom));if (paral_atom) paral_atom=(xmpi_comm_size(comm_atom)>1) 562 nullify(my_atmtab);if (present(mpi_atmtab)) my_atmtab => mpi_atmtab 563 my_comm_atom=xmpi_comm_self;if (present(comm_atom)) my_comm_atom=comm_atom 564 my_atmtab_allocated=.false. 565 566 !Determine in which case we are (parallelism, ...) 567 !No parallelism: a single copy operation 568 paral_case=0;nrhoij_max=nrhoij_in 569 pawrhoij_out => pawrhoij_cpy 570 if (paral_atom) then 571 if (nrhoij_out<nrhoij_in) then ! Parallelism: the copy operation is a scatter 572 call get_my_natom(my_comm_atom,my_nrhoij,nrhoij_in) 573 if (my_nrhoij==nrhoij_out) then 574 call get_my_atmtab(my_comm_atom,my_atmtab,my_atmtab_allocated,paral_atom,nrhoij_in) 575 paral_case=1;nrhoij_max=nrhoij_out 576 pawrhoij_out => pawrhoij_cpy 577 else 578 msg=' nrhoij_out should be equal to my_natom !' 579 LIBPAW_BUG(msg) 580 end if 581 else ! Parallelism: the copy operation is a gather 582 call get_my_natom(my_comm_atom,my_nrhoij,nrhoij_out) 583 if (my_nrhoij==nrhoij_in) then 584 paral_case=2;nrhoij_max=nrhoij_in 585 LIBPAW_DATATYPE_ALLOCATE(pawrhoij_out,(nrhoij_in)) 586 call pawrhoij_nullify(pawrhoij_out) 587 if (nrhoij_in>0) then 588 LIBPAW_ALLOCATE(typat,(nrhoij_in)) 589 LIBPAW_ALLOCATE(nlmn,(nrhoij_in)) 590 do irhoij=1,nrhoij_in 591 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 & use_rhoijres=pawrhoij_cpy(1)%use_rhoijres) 600 LIBPAW_DEALLOCATE(typat) 601 LIBPAW_DEALLOCATE(nlmn) 602 end if 603 else 604 msg=' nrhoij_in should be equal to my_natom!' 605 LIBPAW_BUG(msg) 606 end if 607 end if 608 end if 609 610 !Loop on rhoij components 611 if (nrhoij_max>0) then 612 do irhoij=1,nrhoij_max 613 jrhoij=irhoij;if (paral_case==1) jrhoij=my_atmtab(irhoij) 614 615 lmn2_size_in=pawrhoij_in(jrhoij)%lmn2_size 616 lmn2_size_out=lmn2_size_in 617 cplex_in=pawrhoij_in(jrhoij)%cplex_rhoij 618 cplex_out=cplex_in;if(keep_cplex_)cplex_out=pawrhoij_out(irhoij)%cplex_rhoij 619 qphase_in=pawrhoij_in(jrhoij)%qphase 620 qphase_out=qphase_in;if(keep_qphase_)qphase_out=pawrhoij_out(irhoij)%qphase 621 nspden_in=pawrhoij_in(jrhoij)%nspden 622 nselect=pawrhoij_in(jrhoij)%nrhoijsel 623 nselect_out=pawrhoij_out(irhoij)%nrhoijsel 624 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 & nselect/=nselect_out) 631 cplex=min(cplex_in,cplex_out) 632 qphase=min(qphase_in,qphase_out) 633 634 ! Scalars 635 pawrhoij_out(irhoij)%cplex_rhoij=cplex_out+0 636 pawrhoij_out(irhoij)%qphase=qphase_out+0 637 pawrhoij_out(irhoij)%nspden=nspden_out+0 638 pawrhoij_out(irhoij)%lmn2_size=lmn2_size_out+0 639 pawrhoij_out(irhoij)%lmn_size=pawrhoij_in(jrhoij)%lmn_size+0 640 if(.not.keep_itypat_) pawrhoij_out(irhoij)%itypat =pawrhoij_in(jrhoij)%itypat+0 641 if(.not.keep_nspden_) pawrhoij_out(irhoij)%nsppol =pawrhoij_in(jrhoij)%nsppol+0 642 if(.not.keep_nspden_) pawrhoij_out(irhoij)%nspinor=pawrhoij_in(jrhoij)%nspinor+0 643 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 use_rhoijp=pawrhoij_in(jrhoij)%use_rhoijp 651 if (pawrhoij_out(irhoij)%use_rhoijp/=use_rhoijp) then 652 if (pawrhoij_out(irhoij)%use_rhoijp>0) then 653 LIBPAW_DEALLOCATE(pawrhoij_out(irhoij)%rhoijp) 654 LIBPAW_DEALLOCATE(pawrhoij_out(irhoij)%rhoijselect) 655 end if 656 if (use_rhoijp>0) then 657 LIBPAW_ALLOCATE(pawrhoij_out(irhoij)%rhoijp,(cplex_out*qphase_out*lmn2_size_out,nspden_out)) 658 LIBPAW_ALLOCATE(pawrhoij_out(irhoij)%rhoijselect,(lmn2_size_out)) 659 pawrhoij_out(irhoij)%rhoijselect=0 660 end if 661 end if 662 pawrhoij_out(irhoij)%use_rhoijp=use_rhoijp 663 if (use_rhoijp>0) then 664 if (change_dim) then 665 if(allocated(pawrhoij_out(irhoij)%rhoijp)) then 666 LIBPAW_DEALLOCATE(pawrhoij_out(irhoij)%rhoijp) 667 end if 668 LIBPAW_ALLOCATE(pawrhoij_out(irhoij)%rhoijp,(cplex_out*qphase_out*lmn2_size_out,nspden_out)) 669 end if 670 pawrhoij_out(irhoij)%rhoijp(:,:)=zero 671 if (nspden_out==1) then 672 if (nspden_in==2) then 673 do iphase=1,qphase 674 i_in=(iphase-1)*lmn2_size_in;i_out=(iphase-1)*lmn2_size_out 675 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 & +pawrhoij_in(jrhoij)%rhoijp(i_in+1:i_in+cplex,2)+zero 679 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 do iphase=1,qphase 684 i_in=(iphase-1)*lmn2_size_in;i_out=(iphase-1)*lmn2_size_out 685 do ilmn=1,nselect 686 pawrhoij_out(irhoij)%rhoijp(i_out+1:i_out+cplex,1)= & 687 & pawrhoij_in(jrhoij)%rhoijp(i_in+1:i_in+cplex,1)+zero 688 i_in=i_in+cplex_in;i_out=i_out+cplex_out 689 end do 690 end do 691 end if 692 else if (nspden_out==2) then 693 if (nspden_in==1) then 694 do iphase=1,qphase 695 i_in=(iphase-1)*lmn2_size_in;i_out=(iphase-1)*lmn2_size_out 696 do ilmn=1,nselect 697 pawrhoij_out(irhoij)%rhoijp(i_out+1:i_out+cplex,1)= & 698 & 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 & half*pawrhoij_in(jrhoij)%rhoijp(i_in+1:i_in+cplex,1)+zero 701 i_in=i_in+cplex_in;i_out=i_out+cplex_out 702 end do 703 end do 704 else if (nspden_in==2) then 705 do ispden=1,nspden_out 706 do iphase=1,qphase 707 i_in=(iphase-1)*lmn2_size_in;i_out=(iphase-1)*lmn2_size_out 708 do ilmn=1,nselect 709 pawrhoij_out(irhoij)%rhoijp(i_out+1:i_out+cplex,ispden)= & 710 & pawrhoij_in(jrhoij)%rhoijp(i_in+1:i_in+cplex,ispden)+zero 711 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 do iphase=1,qphase 717 i_in=(iphase-1)*lmn2_size_in;i_out=(iphase-1)*lmn2_size_out 718 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 & +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 & -pawrhoij_in(jrhoij)%rhoijp(i_in+1:i_in+cplex,4))+zero 725 i_in=i_in+cplex_in;i_out=i_out+cplex_out 726 end do 727 end do 728 end if 729 else if (nspden_out==4) then 730 if (nspden_in==1) then 731 do iphase=1,qphase 732 i_in=(iphase-1)*lmn2_size_in;i_out=(iphase-1)*lmn2_size_out 733 do ilmn=1,nselect 734 pawrhoij_out(irhoij)%rhoijp(i_out+1:i_out+cplex,1)= & 735 & pawrhoij_in(jrhoij)%rhoijp(i_in+1:i_in+cplex,1)+zero 736 pawrhoij_out(irhoij)%rhoijp(i_out+1:i_out+cplex,2:4)=zero 737 i_in=i_in+cplex_in;i_out=i_out+cplex_out 738 end do 739 end do 740 else if (nspden_in==2) then 741 do iphase=1,qphase 742 i_in=(iphase-1)*lmn2_size_in;i_out=(iphase-1)*lmn2_size_out 743 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 & +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 & -pawrhoij_in(jrhoij)%rhoijp(i_in+1:i_in+cplex,2)+zero 750 pawrhoij_out(irhoij)%rhoijp(i_out+1:i_out+cplex,2:3)=zero 751 i_in=i_in+cplex_in;i_out=i_out+cplex_out 752 end do 753 end do 754 else ! nspden_in==4 755 do ispden=1,nspden_out 756 do iphase=1,qphase 757 i_in=(iphase-1)*lmn2_size_in;i_out=(iphase-1)*lmn2_size_out 758 do ilmn=1,nselect 759 pawrhoij_out(irhoij)%rhoijp(i_out+1:i_out+cplex,ispden)= & 760 & pawrhoij_in(jrhoij)%rhoijp(i_in+1:i_in+cplex,ispden)+zero 761 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 if (change_dim) then 772 if(allocated(pawrhoij_out(irhoij)%rhoijselect)) then 773 LIBPAW_DEALLOCATE(pawrhoij_out(irhoij)%rhoijselect) 774 end if 775 LIBPAW_ALLOCATE(pawrhoij_out(irhoij)%rhoijselect,(lmn2_size_out)) 776 end if 777 pawrhoij_out(irhoij)%rhoijselect=0 778 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 lmnmix=pawrhoij_in(jrhoij)%lmnmix_sz 783 if (pawrhoij_out(irhoij)%lmnmix_sz/=lmnmix) then 784 if (pawrhoij_out(irhoij)%lmnmix_sz>0) then 785 LIBPAW_DEALLOCATE(pawrhoij_out(irhoij)%kpawmix) 786 end if 787 if (lmnmix>0) then 788 LIBPAW_ALLOCATE(pawrhoij_out(irhoij)%kpawmix,(lmnmix)) 789 pawrhoij_out(irhoij)%kpawmix=0 790 end if 791 pawrhoij_out(irhoij)%lmnmix_sz=lmnmix 792 end if 793 if (lmnmix>0) pawrhoij_out(irhoij)%kpawmix(1:lmnmix)=pawrhoij_in(jrhoij)%kpawmix(1:lmnmix) 794 795 ! Optional pointer: gradients of rhoij 796 ngrhoij=pawrhoij_in(jrhoij)%ngrhoij 797 if (pawrhoij_out(irhoij)%ngrhoij/=ngrhoij) then 798 if (pawrhoij_out(irhoij)%ngrhoij>0) then 799 LIBPAW_DEALLOCATE(pawrhoij_out(irhoij)%grhoij) 800 end if 801 if (ngrhoij>0) then 802 LIBPAW_ALLOCATE(pawrhoij_out(irhoij)%grhoij,(ngrhoij,cplex_out*qphase_out*lmn2_size_out,nspden_out)) 803 end if 804 pawrhoij_out(irhoij)%ngrhoij=ngrhoij 805 end if 806 if (ngrhoij>0) then 807 if (change_dim) then 808 LIBPAW_DEALLOCATE(pawrhoij_out(irhoij)%grhoij) 809 LIBPAW_ALLOCATE(pawrhoij_out(irhoij)%grhoij,(ngrhoij,cplex_out*qphase_out*lmn2_size_out,nspden_out)) 810 end if 811 pawrhoij_out(irhoij)%grhoij(:,:,:)=zero 812 if (nspden_out==1) then 813 if (nspden_in==2) then 814 do iphase=1,qphase 815 i_in=(iphase-1)*lmn2_size_in;i_out=(iphase-1)*lmn2_size_out 816 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 & +pawrhoij_in(jrhoij)%grhoij(1:ngrhoij,i_in+1:i_in+cplex,2)+zero 820 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 do iphase=1,qphase 825 i_in=(iphase-1)*lmn2_size_in;i_out=(iphase-1)*lmn2_size_out 826 do ilmn=1,lmn2_size_out 827 pawrhoij_out(irhoij)%grhoij(1:ngrhoij,i_out+1:i_out+cplex,1)= & 828 & pawrhoij_in(jrhoij)%grhoij(1:ngrhoij,i_in+1:i_in+cplex,1)+zero 829 i_in=i_in+cplex_in;i_out=i_out+cplex_out 830 end do 831 end do 832 end if 833 else if (nspden_out==2) then 834 if (nspden_in==1) then 835 do iphase=1,qphase 836 i_in=(iphase-1)*lmn2_size_in;i_out=(iphase-1)*lmn2_size_out 837 do ilmn=1,lmn2_size_out 838 pawrhoij_out(irhoij)%grhoij(1:ngrhoij,i_out+1:i_out+cplex,1)= & 839 & 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 & half*pawrhoij_in(jrhoij)%grhoij(1:ngrhoij,i_in+1:i_in+cplex,1)+zero 842 i_in=i_in+cplex_in;i_out=i_out+cplex_out 843 end do 844 end do 845 else if (nspden_in==2) then 846 do ispden=1,nspden_out 847 do iphase=1,qphase 848 i_in=(iphase-1)*lmn2_size_in;i_out=(iphase-1)*lmn2_size_out 849 do ilmn=1,lmn2_size_out 850 pawrhoij_out(irhoij)%grhoij(1:ngrhoij,i_out+1:i_out+cplex,ispden)= & 851 & pawrhoij_in(jrhoij)%grhoij(1:ngrhoij,i_in+1:i_in+cplex,ispden)+zero 852 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 do iphase=1,qphase 858 i_in=(iphase-1)*lmn2_size_in;i_out=(iphase-1)*lmn2_size_out 859 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 & +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 & -pawrhoij_in(jrhoij)%grhoij(1:ngrhoij,i_in+1:i_in+cplex,4))+zero 866 i_in=i_in+cplex_in;i_out=i_out+cplex_out 867 end do 868 end do 869 end if 870 else if (nspden_out==4) then 871 if (nspden_in==1) then 872 do iphase=1,qphase 873 i_in=(iphase-1)*lmn2_size_in;i_out=(iphase-1)*lmn2_size_out 874 do ilmn=1,lmn2_size_out 875 pawrhoij_out(irhoij)%grhoij(1:ngrhoij,i_out+1:i_out+cplex,1)= & 876 & pawrhoij_in(jrhoij)%grhoij(1:ngrhoij,i_in+1:i_in+cplex,1)+zero 877 pawrhoij_out(irhoij)%grhoij(1:ngrhoij,i_out+1:i_out+cplex,2:4)=zero 878 i_in=i_in+cplex_in;i_out=i_out+cplex_out 879 end do 880 end do 881 else if (nspden_in==2) then 882 do iphase=1,qphase 883 i_in=(iphase-1)*lmn2_size_in;i_out=(iphase-1)*lmn2_size_out 884 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 & +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 & -pawrhoij_in(jrhoij)%grhoij(1:ngrhoij,i_in+1:i_in+cplex,2)+zero 891 pawrhoij_out(irhoij)%grhoij(1:ngrhoij,i_out+1:i_out+cplex,2:3)=zero 892 i_in=i_in+cplex_in;i_out=i_out+cplex_out 893 end do 894 end do 895 else ! nspden_in==4 896 do ispden=1,nspden_out 897 do iphase=1,qphase 898 i_in=(iphase-1)*lmn2_size_in;i_out=(iphase-1)*lmn2_size_out 899 do ilmn=1,lmn2_size_out 900 pawrhoij_out(irhoij)%grhoij(1:ngrhoij,i_out+1:i_out+cplex,ispden)= & 901 & pawrhoij_in(jrhoij)%grhoij(1:ngrhoij,i_in+1:i_in+cplex,ispden)+zero 902 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 use_rhoijres=pawrhoij_in(jrhoij)%use_rhoijres 912 if (pawrhoij_out(irhoij)%use_rhoijres/=use_rhoijres) then 913 if (pawrhoij_out(irhoij)%use_rhoijres>0) then 914 LIBPAW_DEALLOCATE(pawrhoij_out(irhoij)%rhoijres) 915 end if 916 if (use_rhoijres>0) then 917 LIBPAW_ALLOCATE(pawrhoij_out(irhoij)%rhoijres,(cplex_out*qphase_out*lmn2_size_out,nspden_out)) 918 end if 919 pawrhoij_out(irhoij)%use_rhoijres=use_rhoijres 920 end if 921 if (use_rhoijres>0) then 922 if (change_dim) then 923 LIBPAW_DEALLOCATE(pawrhoij_out(irhoij)%rhoijres) 924 LIBPAW_ALLOCATE(pawrhoij_out(irhoij)%rhoijres,(cplex_out*qphase_out*lmn2_size_out,nspden_out)) 925 end if 926 pawrhoij_out(irhoij)%rhoijres(:,:)=zero 927 if (nspden_out==1) then 928 if (nspden_in==2) then 929 do iphase=1,qphase 930 i_in=(iphase-1)*lmn2_size_in;i_out=(iphase-1)*lmn2_size_out 931 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 & +pawrhoij_in(jrhoij)%rhoijres(i_in+1:i_in+cplex,2)+zero 935 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 do iphase=1,qphase 940 i_in=(iphase-1)*lmn2_size_in;i_out=(iphase-1)*lmn2_size_out 941 do ilmn=1,lmn2_size_out 942 pawrhoij_out(irhoij)%rhoijres(i_out+1:i_out+cplex,1)= & 943 & pawrhoij_in(jrhoij)%rhoijres(i_in+1:i_in+cplex,1)+zero 944 i_in=i_in+cplex_in;i_out=i_out+cplex_out 945 end do 946 end do 947 end if 948 else if (nspden_out==2) then 949 if (nspden_in==1) then 950 do iphase=1,qphase 951 i_in=(iphase-1)*lmn2_size_in;i_out=(iphase-1)*lmn2_size_out 952 do ilmn=1,lmn2_size_out 953 pawrhoij_out(irhoij)%rhoijres(i_out+1:i_out+cplex,1)= & 954 & 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 & half*pawrhoij_in(jrhoij)%rhoijres(i_in+1:i_in+cplex,1)+zero 957 i_in=i_in+cplex_in;i_out=i_out+cplex_out 958 end do 959 end do 960 else if (nspden_in==2) then 961 do ispden=1,nspden_out 962 do iphase=1,qphase 963 i_in=(iphase-1)*lmn2_size_in;i_out=(iphase-1)*lmn2_size_out 964 do ilmn=1,lmn2_size_out 965 pawrhoij_out(irhoij)%rhoijres(i_out+1:i_out+cplex,ispden)= & 966 & pawrhoij_in(jrhoij)%rhoijres(i_in+1:i_in+cplex,ispden)+zero 967 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 do iphase=1,qphase 973 i_in=(iphase-1)*lmn2_size_in;i_out=(iphase-1)*lmn2_size_out 974 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 & +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 & -pawrhoij_in(jrhoij)%rhoijres(i_in+1:i_in+cplex,4))+zero 981 i_in=i_in+cplex_in;i_out=i_out+cplex_out 982 end do 983 end do 984 end if 985 else if (nspden_out==4) then 986 if (nspden_in==1) then 987 do iphase=1,qphase 988 i_in=(iphase-1)*lmn2_size_in;i_out=(iphase-1)*lmn2_size_out 989 do ilmn=1,lmn2_size_out 990 pawrhoij_out(irhoij)%rhoijres(i_out+1:i_out+cplex,1)= & 991 & pawrhoij_in(jrhoij)%rhoijres(i_in+1:i_in+cplex,1)+zero 992 pawrhoij_out(irhoij)%rhoijres(i_out+1:i_out+cplex,2:4)=zero 993 i_in=i_in+cplex_in;i_out=i_out+cplex_out 994 end do 995 end do 996 else if (nspden_in==2) then 997 do iphase=1,qphase 998 i_in=(iphase-1)*lmn2_size_in;i_out=(iphase-1)*lmn2_size_out 999 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 & +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 & -pawrhoij_in(jrhoij)%rhoijres(i_in+1:i_in+cplex,2)+zero 1006 pawrhoij_out(irhoij)%rhoijres(i_out+1:i_out+cplex,2:3)=zero 1007 i_in=i_in+cplex_in;i_out=i_out+cplex_out 1008 end do 1009 end do 1010 else ! nspden_in==4 1011 do ispden=1,nspden_out 1012 do iphase=1,qphase 1013 i_in=(iphase-1)*lmn2_size_in;i_out=(iphase-1)*lmn2_size_out 1014 do ilmn=1,lmn2_size_out 1015 pawrhoij_out(irhoij)%rhoijres(i_out+1:i_out+cplex,ispden)= & 1016 & pawrhoij_in(jrhoij)%rhoijres(i_in+1:i_in+cplex,ispden)+zero 1017 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 use_rhoij_=pawrhoij_in(jrhoij)%use_rhoij_ 1027 if (pawrhoij_out(irhoij)%use_rhoij_/=use_rhoij_) then 1028 if (pawrhoij_out(irhoij)%use_rhoij_>0) then 1029 LIBPAW_DEALLOCATE(pawrhoij_out(irhoij)%rhoij_) 1030 end if 1031 if (use_rhoij_>0) then 1032 LIBPAW_ALLOCATE(pawrhoij_out(irhoij)%rhoij_,(cplex_out*qphase_out*lmn2_size_out,nspden_out)) 1033 end if 1034 pawrhoij_out(irhoij)%use_rhoij_=use_rhoij_ 1035 end if 1036 if (use_rhoij_>0) then 1037 if (change_dim) then 1038 if(allocated(pawrhoij_out(irhoij)%rhoij_)) then 1039 LIBPAW_DEALLOCATE(pawrhoij_out(irhoij)%rhoij_) 1040 end if 1041 LIBPAW_ALLOCATE(pawrhoij_out(irhoij)%rhoij_,(cplex_out*qphase_out*lmn2_size_out,nspden_out)) 1042 end if 1043 pawrhoij_out(irhoij)%rhoij_(:,:)=zero 1044 if (nspden_out==1) then 1045 if (nspden_in==2) then 1046 do iphase=1,qphase 1047 i_in=(iphase-1)*lmn2_size_in;i_out=(iphase-1)*lmn2_size_out 1048 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 & +pawrhoij_in(jrhoij)%rhoij_(i_in+1:i_in+cplex,2)+zero 1052 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 do iphase=1,qphase 1057 i_in=(iphase-1)*lmn2_size_in;i_out=(iphase-1)*lmn2_size_out 1058 do ilmn=1,lmn2_size_out 1059 pawrhoij_out(irhoij)%rhoij_(i_out+1:i_out+cplex,1)= & 1060 & pawrhoij_in(jrhoij)%rhoij_(i_in+1:i_in+cplex,1)+zero 1061 i_in=i_in+cplex_in;i_out=i_out+cplex_out 1062 end do 1063 end do 1064 end if 1065 else if (nspden_out==2) then 1066 if (nspden_in==1) then 1067 do iphase=1,qphase 1068 i_in=(iphase-1)*lmn2_size_in;i_out=(iphase-1)*lmn2_size_out 1069 do ilmn=1,lmn2_size_out 1070 pawrhoij_out(irhoij)%rhoij_(i_out+1:i_out+cplex,1)= & 1071 & 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 & half*pawrhoij_in(irhoij)%rhoij_(i_in+1:i_in+cplex,1)+zero 1074 i_in=i_in+cplex_in;i_out=i_out+cplex_out 1075 end do 1076 end do 1077 else if (nspden_in==2) then 1078 do ispden=1,nspden_out 1079 do iphase=1,qphase 1080 i_in=(iphase-1)*lmn2_size_in;i_out=(iphase-1)*lmn2_size_out 1081 do ilmn=1,lmn2_size_out 1082 pawrhoij_out(irhoij)%rhoij_(i_out+1:i_out+cplex,ispden)= & 1083 & pawrhoij_in(jrhoij)%rhoij_(i_in+1:i_in+cplex,ispden)+zero 1084 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 do iphase=1,qphase 1090 i_in=(iphase-1)*lmn2_size_in;i_out=(iphase-1)*lmn2_size_out 1091 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 & +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 & -pawrhoij_in(jrhoij)%rhoij_(i_in+1:i_in+cplex,4))+zero 1098 i_in=i_in+cplex_in;i_out=i_out+cplex_out 1099 end do 1100 end do 1101 end if 1102 else if (nspden_out==4) then 1103 if (nspden_in==1) then 1104 do iphase=1,qphase 1105 i_in=(iphase-1)*lmn2_size_in;i_out=(iphase-1)*lmn2_size_out 1106 do ilmn=1,lmn2_size_out 1107 pawrhoij_out(irhoij)%rhoij_(i_out+1:i_out+cplex,1)= & 1108 & pawrhoij_in(jrhoij)%rhoij_(i_in+1:i_in+cplex,1)+zero 1109 pawrhoij_out(irhoij)%rhoij_(i_out+1:i_out+cplex,2:4)=zero 1110 i_in=i_in+cplex_in;i_out=i_out+cplex_out 1111 end do 1112 end do 1113 else if (nspden_in==2) then 1114 do iphase=1,qphase 1115 i_in=(iphase-1)*lmn2_size_in;i_out=(iphase-1)*lmn2_size_out 1116 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 & +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 & -pawrhoij_in(jrhoij)%rhoij_(i_in+1:i_in+cplex,2)+zero 1123 pawrhoij_out(irhoij)%rhoij_(i_out+1:i_out+cplex,2:3)=zero 1124 i_in=i_in+cplex_in;i_out=i_out+cplex_out 1125 end do 1126 end do 1127 else ! nspden_in==4 1128 do ispden=1,nspden_out 1129 do iphase=1,qphase 1130 i_in=(iphase-1)*lmn2_size_in;i_out=(iphase-1)*lmn2_size_out 1131 do ilmn=1,lmn2_size_out 1132 pawrhoij_out(irhoij)%rhoij_(i_out+1:i_out+cplex,ispden)= & 1133 & pawrhoij_in(jrhoij)%rhoij_(i_in+1:i_in+cplex,ispden)+zero 1134 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 if (paral_case==2) then 1147 call pawrhoij_free(pawrhoij_cpy) 1148 call pawrhoij_gather(pawrhoij_out,pawrhoij_cpy,-1,my_comm_atom) 1149 call pawrhoij_free(pawrhoij_out) 1150 LIBPAW_DATATYPE_DEALLOCATE(pawrhoij_out) 1151 1152 ! Sequential case: fill missing elements 1153 else if (paral_case==0) then 1154 if (nrhoij_in<nrhoij_out) then 1155 do irhoij=nrhoij_in+1,nrhoij_out 1156 pawrhoij_cpy(irhoij)%nrhoijsel=0 1157 if (pawrhoij_cpy(irhoij)%use_rhoijp>0) then 1158 pawrhoij_cpy(irhoij)%rhoijselect=0 1159 pawrhoij_cpy(irhoij)%rhoijp=zero 1160 end if 1161 if (pawrhoij_cpy(irhoij)%lmnmix_sz>0) pawrhoij_cpy(irhoij)%kpawmix=0 1162 if (pawrhoij_cpy(irhoij)%ngrhoij>0) pawrhoij_cpy(irhoij)%grhoij=zero 1163 if (pawrhoij_cpy(irhoij)%use_rhoij_>0) pawrhoij_cpy(irhoij)%rhoij_=zero 1164 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 call free_my_atmtab(my_atmtab,my_atmtab_allocated) 1171 1172 end subroutine pawrhoij_copy
m_pawrhoij/pawrhoij_filter [ Functions ]
[ Top ] [ m_pawrhoij ] [ Functions ]
NAME
pawrhoij_filter
FUNCTION
Filter a "rhoij" array (PAW on-site occupancies), i.e. select only the non-zero elements. INPUT cplex=1 if Rhoij is real, 2 if Rhoij is complex qphase=2 if rhoij has a exp(iqR) phase, 1 if not lmn2_size=dimension of rhoij=number of (i,j) pairs with i<=j nspden=number of rhoij spin components [rhoij_input(cplex*qphase*lmn2_size,nspden)]= -- optional argument-- input values for rhoij (including zero values) If this argument is not provided, that the input values from rhoij()
OUTPUT
nselect=number of non-zero values of rhoij
rhoijselect(lmn2_size)=contains the indices of the selected (i,j) pairs
rhoijselect(nselect+1:lmn2_size)=0
SIDE EFFECTS
rhoij(cplex*qphase*lmn2_size,nspden)=
* Input: the array is filled with all rhoij values (only if rhoij_input is not present)
* Output: the nselect first elements contain the non-zero rhoij values,
next value are irrelevant
SOURCE
3083 subroutine pawrhoij_filter(rhoij,rhoijselect,nselect,cplex,qphase,lmn2_size,nspden, & 3084 & 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 real(dp),pointer :: rhoij_in(:,:) 3101 3102 ! ************************************************************************* 3103 3104 nselect=0 3105 rhoijselect(:)=0 3106 3107 if (present(rhoij_input)) then 3108 rhoij_in => rhoij_input 3109 else 3110 rhoij_in => rhoij 3111 end if 3112 3113 !Treat each cplex/qphase case separately 3114 3115 if (cplex==1) then 3116 3117 if (qphase==1) then 3118 3119 do klmn=1,lmn2_size 3120 if (any(abs(rhoij_in(klmn,:))>tol_rhoij)) then 3121 nselect=nselect+1 3122 rhoijselect(nselect)=klmn 3123 do isp=1,nspden 3124 rhoij(nselect,isp)=rhoij_in(klmn,isp) 3125 end do 3126 end if 3127 end do 3128 3129 else if (qphase==2) then 3130 3131 do klmn=1,lmn2_size 3132 klmn2=klmn+lmn2_size 3133 if (any(abs(rhoij_in(klmn,:))>tol_rhoij).or. & 3134 & any(abs(rhoij_in(klmn2,:))>tol_rhoij)) then 3135 nselect=nselect+1 ; nsel2=nselect+lmn2_size 3136 rhoijselect(nselect)=klmn 3137 do isp=1,nspden 3138 rhoij(nselect,isp)=rhoij_in(klmn ,isp) 3139 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 if (qphase==1) then 3149 do klmn=1,lmn2_size 3150 klmn1=2*klmn 3151 if (any(abs(rhoij_in(klmn1-1:klmn1,:))>tol_rhoij)) then 3152 nselect=nselect+1 ; nsel1=2*nselect 3153 rhoijselect(nselect)=klmn 3154 do isp=1,nspden 3155 rhoij(nsel1-1,isp)=rhoij_in(klmn1-1,isp) 3156 rhoij(nsel1 ,isp)=rhoij_in(klmn1 ,isp) 3157 end do 3158 end if 3159 end do 3160 3161 else if (qphase==2) then 3162 3163 do klmn=1,lmn2_size 3164 klmn1=2*klmn ; klmn2=klmn1+lmn2_size 3165 if (any(abs(rhoij_in(klmn1-1:klmn1,:))>tol_rhoij).or. & 3166 & any(abs(rhoij_in(klmn2-1:klmn2,:))>tol_rhoij)) then 3167 nselect=nselect+1 ; nsel1=2*nselect ; nsel2=nsel1+lmn2_size 3168 rhoijselect(nselect)=klmn 3169 do isp=1,nspden 3170 rhoij(nsel1-1,isp)=rhoij_in(klmn1-1,isp) 3171 rhoij(nsel1 ,isp)=rhoij_in(klmn1 ,isp) 3172 rhoij(nsel2-1,isp)=rhoij_in(klmn2-1,isp) 3173 rhoij(nsel2 ,isp)=rhoij_in(klmn2 ,isp) 3174 end do 3175 end if 3176 end do 3177 3178 end if 3179 endif 3180 3181 end subroutine pawrhoij_filter
m_pawrhoij/pawrhoij_free [ Functions ]
[ Top ] [ m_pawrhoij ] [ Functions ]
NAME
pawrhoij_free
FUNCTION
Destroy a pawrhoij datastructure
SIDE EFFECTS
pawrhoij(:)<type(pawrhoij_type)>= rhoij datastructure
SOURCE
386 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 nrhoij=size(pawrhoij) 399 400 if (nrhoij>0) then 401 do irhoij=1,nrhoij 402 pawrhoij(irhoij)%cplex_rhoij=1 403 pawrhoij(irhoij)%qphase=1 404 pawrhoij(irhoij)%nrhoijsel=0 405 pawrhoij(irhoij)%ngrhoij=0 406 pawrhoij(irhoij)%lmnmix_sz=0 407 pawrhoij(irhoij)%use_rhoij_=0 408 pawrhoij(irhoij)%use_rhoijp=0 409 pawrhoij(irhoij)%use_rhoijres=0 410 if (allocated(pawrhoij(irhoij)%rhoijp)) then 411 LIBPAW_DEALLOCATE(pawrhoij(irhoij)%rhoijp) 412 end if 413 if (allocated(pawrhoij(irhoij)%rhoijselect)) then 414 LIBPAW_DEALLOCATE(pawrhoij(irhoij)%rhoijselect) 415 end if 416 if (allocated(pawrhoij(irhoij)%grhoij)) then 417 LIBPAW_DEALLOCATE(pawrhoij(irhoij)%grhoij) 418 end if 419 if (allocated(pawrhoij(irhoij)%kpawmix)) then 420 LIBPAW_DEALLOCATE(pawrhoij(irhoij)%kpawmix) 421 end if 422 if (allocated(pawrhoij(irhoij)%rhoij_)) then 423 LIBPAW_DEALLOCATE(pawrhoij(irhoij)%rhoij_) 424 end if 425 if (allocated(pawrhoij(irhoij)%rhoijres)) then 426 LIBPAW_DEALLOCATE(pawrhoij(irhoij)%rhoijres) 427 end if 428 end do 429 end if 430 431 end subroutine pawrhoij_free
m_pawrhoij/pawrhoij_free_unpacked [ Functions ]
[ Top ] [ m_pawrhoij ] [ Functions ]
NAME
pawrhoij_free_unpacked
FUNCTION
Destroy field of rhoij datastructure for unpacked values (pawrhoij%rhoij_ array)
SIDE EFFECTS
rhoij(:) <pawrhoij_type)>= input/output datastructure * In output the rhoij_ array is deallocated
SOURCE
2848 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 nrhoij=SIZE(rhoij);if (nrhoij==0) return 2861 2862 do iat=1,nrhoij 2863 2864 if (allocated(rhoij(iat)%rhoij_)) then 2865 LIBPAW_DEALLOCATE(rhoij(iat)%rhoij_) 2866 end if 2867 rhoij(iat)%use_rhoij_=0 2868 2869 end do 2870 2871 end subroutine pawrhoij_free_unpacked
m_pawrhoij/pawrhoij_gather [ Functions ]
[ Top ] [ m_pawrhoij ] [ Functions ]
NAME
pawrhoij_gather
FUNCTION
(All)Gather pawrhoij datastructures
INPUTS
master=master communicator receiving data ; if -1 do a ALLGATHER comm_atom= communicator pawrhoij_in(:)<type(pawrhoij_type)>= input rhoij datastructures on every process with_grhoij : optional argument (logical, default=.TRUE.) TRUE if pawrhoij%grhoij field is included in the gather operation with_lmnmix : optional argument (logical, default=.TRUE.) TRUE if pawrhoij%lmnmix field is included in the gather operation with_rhoijp : optional argument (logical, default=.TRUE.) TRUE if pawrhoij%rhoijp and pawrhoij%rhoijselect fields are included in the gather operation with_rhoijres : optional argument (logical, default=.TRUE.) TRUE if pawrhoij%rhoijres field is included in the gather operation with_rhoij_ : optional argument (logical, default=.TRUE.) TRUE if pawrhoij%rhoij_ field is included in the gather operation
OUTPUT
pawrhoij_gathered(:)<type(pawrhoij_type)>= output rhoij datastructure
NOTES
The gathered structure are ordered like in sequential mode.
SOURCE
1209 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 integer,allocatable :: buf_int(:),buf_int_all(:) 1231 integer,allocatable :: count_dp(:),count_int(:),count_tot(:),displ_dp(:),displ_int(:) 1232 integer, pointer :: my_atmtab(:) 1233 real(dp),allocatable :: buf_dp(:),buf_dp_all(:) 1234 1235 ! ************************************************************************* 1236 1237 nrhoij_in=size(pawrhoij_in);nrhoij_out=size(pawrhoij_gathered) 1238 1239 nproc_atom=xmpi_comm_size(comm_atom) 1240 me_atom=xmpi_comm_rank(comm_atom) 1241 1242 if (nproc_atom==1) then 1243 if (master==-1.or.me_atom==master) then 1244 call pawrhoij_copy(pawrhoij_in,pawrhoij_gathered,.false.,.false.,.false.) 1245 end if 1246 return 1247 end if 1248 1249 !Test on sizes 1250 nrhoij_in_sum=nrhoij_in 1251 call xmpi_sum(nrhoij_in_sum,comm_atom,ierr) 1252 if (master==-1) then 1253 if (nrhoij_out/=nrhoij_in_sum) then 1254 msg='Wrong sizes sum[nrhoij_ij]/=nrhoij_out !' 1255 LIBPAW_BUG(msg) 1256 end if 1257 else 1258 if (me_atom==master.and.nrhoij_out/=nrhoij_in_sum) then 1259 msg='(2) pawrhoij_gathered wrongly allocated !' 1260 LIBPAW_BUG(msg) 1261 end if 1262 end if 1263 1264 !Optional arguments 1265 with_grhoij_ =.true.;if (present(with_grhoij)) with_grhoij_ =with_grhoij 1266 with_lmnmix_ =.true.;if (present(with_lmnmix)) with_lmnmix_ =with_lmnmix 1267 with_rhoijp_ =.true.;if (present(with_rhoijp)) with_rhoijp_ =with_rhoijp 1268 with_rhoijres_=.true.;if (present(with_rhoijres))with_rhoijres_=with_rhoijres 1269 with_rhoij__ =.true.;if (present(with_rhoij_)) with_rhoij__ =with_rhoij_ 1270 1271 !Retrieve table of atoms 1272 paral_atom=.true.;nullify(my_atmtab) 1273 call get_my_atmtab(comm_atom,my_atmtab,my_atmtab_allocated,paral_atom,nrhoij_in_sum, & 1274 & my_natom_ref=nrhoij_in) 1275 1276 !Compute sizes of buffers 1277 buf_int_size=0;buf_dp_size=0 1278 nselect=0;lmnmix=0;ngrhoij=0;rhoij_size2=0 1279 use_rhoijp=0;use_rhoijres=0;use_rhoij_=0 1280 do irhoij=1,nrhoij_in 1281 cplex =pawrhoij_in(irhoij)%cplex_rhoij 1282 qphase =pawrhoij_in(irhoij)%qphase 1283 lmn2_size=pawrhoij_in(irhoij)%lmn2_size 1284 nspden =pawrhoij_in(irhoij)%nspden 1285 if (with_lmnmix_) lmnmix=pawrhoij_in(irhoij)%lmnmix_sz 1286 if (with_grhoij_) ngrhoij=pawrhoij_in(irhoij)%ngrhoij 1287 if (with_rhoijp_) use_rhoijp=pawrhoij_in(irhoij)%use_rhoijp 1288 if (with_rhoijres_)use_rhoijres=pawrhoij_in(irhoij)%use_rhoijres 1289 if (with_rhoij__) use_rhoij_=pawrhoij_in(irhoij)%use_rhoij_ 1290 buf_int_size=buf_int_size+16 1291 if (use_rhoijp>0) then 1292 nselect=pawrhoij_in(irhoij)%nrhoijsel 1293 buf_int_size=buf_int_size+nselect 1294 buf_dp_size=buf_dp_size + cplex*qphase*nselect*nspden 1295 end if 1296 if (lmnmix>0) buf_int_size=buf_int_size+lmnmix 1297 if (ngrhoij>0) buf_dp_size=buf_dp_size + cplex*qphase*lmn2_size*nspden*ngrhoij 1298 if (use_rhoijres>0) buf_dp_size=buf_dp_size + cplex*qphase*lmn2_size*nspden 1299 if (use_rhoij_>0) then 1300 rhoij_size2=size(pawrhoij_in(irhoij)%rhoij_,dim=2) 1301 buf_dp_size=buf_dp_size + cplex*qphase*lmn2_size*rhoij_size2 1302 end if 1303 end do 1304 1305 !Fill input buffers 1306 LIBPAW_ALLOCATE(buf_int,(buf_int_size)) 1307 LIBPAW_ALLOCATE(buf_dp ,(buf_dp_size)) 1308 indx_int=1;indx_dp =1 1309 lmnmix=0;ngrhoij=0;nselect=0;rhoij_size2=0 1310 use_rhoijp=0;use_rhoijres=0;use_rhoij_=0 1311 do irhoij=1,nrhoij_in 1312 cplex =pawrhoij_in(irhoij)%cplex_rhoij 1313 qphase =pawrhoij_in(irhoij)%qphase 1314 lmn2_size=pawrhoij_in(irhoij)%lmn2_size 1315 nspden =pawrhoij_in(irhoij)%nspden 1316 if (with_lmnmix_) lmnmix=pawrhoij_in(irhoij)%lmnmix_sz 1317 if (with_grhoij_) ngrhoij=pawrhoij_in(irhoij)%ngrhoij 1318 if (with_rhoijp_) use_rhoijp=pawrhoij_in(irhoij)%use_rhoijp 1319 if (use_rhoijp > 0) nselect=pawrhoij_in(irhoij)%nrhoijsel 1320 if (with_rhoijres_)use_rhoijres=pawrhoij_in(irhoij)%use_rhoijres 1321 if (with_rhoij__) use_rhoij_ =pawrhoij_in(irhoij)%use_rhoij_ 1322 if (use_rhoij_> 0) rhoij_size2 =size(pawrhoij_in(irhoij)%rhoij_,dim=2) 1323 buf_int(indx_int)=my_atmtab(irhoij) ;indx_int=indx_int+1 1324 buf_int(indx_int)=cplex ;indx_int=indx_int+1 1325 buf_int(indx_int)=qphase ;indx_int=indx_int+1 1326 buf_int(indx_int)=lmn2_size ;indx_int=indx_int+1 1327 buf_int(indx_int)=nspden ;indx_int=indx_int+1 1328 buf_int(indx_int)=nselect ;indx_int=indx_int+1 1329 buf_int(indx_int)=lmnmix ;indx_int=indx_int+1 1330 buf_int(indx_int)=ngrhoij ;indx_int=indx_int+1 1331 buf_int(indx_int)=use_rhoijp ;indx_int=indx_int+1 1332 buf_int(indx_int)=use_rhoijres ;indx_int=indx_int+1 1333 buf_int(indx_int)=use_rhoij_ ;indx_int=indx_int+1 1334 buf_int(indx_int)=rhoij_size2 ;indx_int=indx_int+1 1335 buf_int(indx_int)=pawrhoij_in(irhoij)%itypat ;indx_int=indx_int+1 1336 buf_int(indx_int)=pawrhoij_in(irhoij)%lmn_size ;indx_int=indx_int+1 1337 buf_int(indx_int)=pawrhoij_in(irhoij)%nsppol ;indx_int=indx_int+1 1338 buf_int(indx_int)=pawrhoij_in(irhoij)%nspinor ;indx_int=indx_int+1 1339 if (use_rhoijp>0) then 1340 buf_int(indx_int:indx_int+nselect-1)=pawrhoij_in(irhoij)%rhoijselect(1:nselect) 1341 indx_int=indx_int+nselect 1342 do isp=1,nspden 1343 do ii=1,qphase 1344 jj=(ii-1)*cplex*lmn2_size 1345 buf_dp(indx_dp:indx_dp+cplex*nselect-1)= & 1346 & pawrhoij_in(irhoij)%rhoijp(jj+1:jj+cplex*nselect,isp) 1347 indx_dp=indx_dp+cplex*nselect 1348 end do 1349 end do 1350 end if 1351 if (lmnmix>0) then 1352 buf_int(indx_int:indx_int+lmnmix-1)=pawrhoij_in(irhoij)%kpawmix(1:lmnmix) 1353 indx_int=indx_int+lmnmix 1354 end if 1355 if (ngrhoij>0) then 1356 do isp=1,nspden 1357 do ii=1,cplex*qphase*lmn2_size 1358 buf_dp(indx_dp:indx_dp+ngrhoij-1)=pawrhoij_in(irhoij)%grhoij(1:ngrhoij,ii,isp) 1359 indx_dp=indx_dp+ngrhoij 1360 end do 1361 end do 1362 end if 1363 if (use_rhoijres>0) then 1364 do isp=1,nspden 1365 buf_dp(indx_dp:indx_dp+cplex*qphase*lmn2_size-1)= & 1366 pawrhoij_in(irhoij)%rhoijres(1:cplex*qphase*lmn2_size,isp) 1367 indx_dp=indx_dp+cplex*qphase*lmn2_size 1368 end do 1369 end if 1370 if (use_rhoij_>0) then 1371 do isp=1,rhoij_size2 1372 buf_dp(indx_dp:indx_dp+cplex*qphase*lmn2_size-1)= & 1373 & pawrhoij_in(irhoij)%rhoij_(1:cplex*qphase*lmn2_size,isp) 1374 indx_dp=indx_dp+cplex*qphase*lmn2_size 1375 end do 1376 end if 1377 end do 1378 1379 !Check 1380 if ((indx_int-1/=buf_int_size).or.(indx_dp-1/=buf_dp_size)) then 1381 write(msg,*) 'Wrong buffer sizes: buf_int_size=',buf_int_size,' buf_dp_size=',buf_dp_size 1382 LIBPAW_BUG(msg) 1383 end if 1384 1385 !Communicate (1 gather for integers, 1 gather for reals) 1386 LIBPAW_ALLOCATE(count_int,(nproc_atom)) 1387 LIBPAW_ALLOCATE(displ_int,(nproc_atom)) 1388 LIBPAW_ALLOCATE(count_dp ,(nproc_atom)) 1389 LIBPAW_ALLOCATE(displ_dp ,(nproc_atom)) 1390 LIBPAW_ALLOCATE(count_tot,(2*nproc_atom)) 1391 bufsz(1)=buf_int_size; bufsz(2)=buf_dp_size 1392 call xmpi_allgather(bufsz,2,count_tot,comm_atom,ierr) 1393 do ii=1,nproc_atom 1394 count_int(ii)=count_tot(2*ii-1) 1395 count_dp (ii)=count_tot(2*ii) 1396 end do 1397 displ_int(1)=0;displ_dp(1)=0 1398 do ii=2,nproc_atom 1399 displ_int(ii)=displ_int(ii-1)+count_int(ii-1) 1400 displ_dp (ii)=displ_dp (ii-1)+count_dp (ii-1) 1401 end do 1402 buf_int_size_all=sum(count_int) 1403 buf_dp_size_all =sum(count_dp) 1404 LIBPAW_DEALLOCATE(count_tot) 1405 if (master==-1.or.me_atom==master) then 1406 LIBPAW_ALLOCATE(buf_int_all,(buf_int_size_all)) 1407 LIBPAW_ALLOCATE(buf_dp_all ,(buf_dp_size_all)) 1408 else 1409 LIBPAW_ALLOCATE(buf_int_all,(0)) 1410 LIBPAW_ALLOCATE(buf_dp_all ,(0)) 1411 end if 1412 if (master==-1) then 1413 call xmpi_allgatherv(buf_int,buf_int_size,buf_int_all,count_int,displ_int,comm_atom,ierr) 1414 call xmpi_allgatherv(buf_dp ,buf_dp_size ,buf_dp_all ,count_dp ,displ_dp ,comm_atom,ierr) 1415 else 1416 call xmpi_gatherv(buf_int,buf_int_size,buf_int_all,count_int,displ_int,master,comm_atom,ierr) 1417 call xmpi_gatherv(buf_dp ,buf_dp_size ,buf_dp_all ,count_dp ,displ_dp ,master,comm_atom,ierr) 1418 end if 1419 LIBPAW_DEALLOCATE(count_int) 1420 LIBPAW_DEALLOCATE(displ_int) 1421 LIBPAW_DEALLOCATE(count_dp) 1422 LIBPAW_DEALLOCATE(displ_dp) 1423 1424 !Retrieve data from output buffer 1425 if (master==-1.or.me_atom==master) then 1426 indx_int=1;indx_dp=1 1427 call pawrhoij_free(pawrhoij_gathered) 1428 do irhoij=1,nrhoij_out 1429 jrhoij =buf_int_all(indx_int) ;indx_int=indx_int+1 1430 cplex =buf_int_all(indx_int) ;indx_int=indx_int+1 1431 qphase =buf_int_all(indx_int) ;indx_int=indx_int+1 1432 lmn2_size =buf_int_all(indx_int) ;indx_int=indx_int+1 1433 nspden =buf_int_all(indx_int) ;indx_int=indx_int+1 1434 nselect =buf_int_all(indx_int) ;indx_int=indx_int+1 1435 lmnmix =buf_int_all(indx_int) ;indx_int=indx_int+1 1436 ngrhoij =buf_int_all(indx_int) ;indx_int=indx_int+1 1437 use_rhoijp =buf_int_all(indx_int) ;indx_int=indx_int+1 1438 use_rhoijres=buf_int_all(indx_int) ;indx_int=indx_int+1 1439 use_rhoij_ =buf_int_all(indx_int) ;indx_int=indx_int+1 1440 rhoij_size2 =buf_int_all(indx_int) ;indx_int=indx_int+1 1441 pawrhoij_gathered(jrhoij)%itypat=buf_int_all(indx_int) ;indx_int=indx_int+1 1442 pawrhoij_gathered(jrhoij)%lmn_size=buf_int_all(indx_int) ;indx_int=indx_int+1 1443 pawrhoij_gathered(jrhoij)%nsppol=buf_int_all(indx_int) ;indx_int=indx_int+1 1444 pawrhoij_gathered(jrhoij)%nspinor=buf_int_all(indx_int) ;indx_int=indx_int+1 1445 pawrhoij_gathered(jrhoij)%cplex_rhoij=cplex 1446 pawrhoij_gathered(jrhoij)%qphase=qphase 1447 pawrhoij_gathered(jrhoij)%lmn2_size=lmn2_size 1448 pawrhoij_gathered(jrhoij)%nspden=nspden 1449 pawrhoij_gathered(jrhoij)%nrhoijsel=nselect 1450 pawrhoij_gathered(jrhoij)%lmnmix_sz=lmnmix 1451 pawrhoij_gathered(jrhoij)%ngrhoij=ngrhoij 1452 pawrhoij_gathered(jrhoij)%use_rhoijp=use_rhoijp 1453 pawrhoij_gathered(jrhoij)%use_rhoijres=use_rhoijres 1454 pawrhoij_gathered(jrhoij)%use_rhoij_=use_rhoij_ 1455 if (use_rhoijp>0) then 1456 LIBPAW_ALLOCATE(pawrhoij_gathered(jrhoij)%rhoijselect,(lmn2_size)) 1457 pawrhoij_gathered(jrhoij)%rhoijselect(1:nselect)=buf_int_all(indx_int:indx_int+nselect-1) 1458 if (nselect < lmn2_size )pawrhoij_gathered(jrhoij)%rhoijselect(nselect+1:lmn2_size)=0 1459 indx_int=indx_int+nselect 1460 LIBPAW_ALLOCATE(pawrhoij_gathered(jrhoij)%rhoijp,(qphase*cplex*lmn2_size,nspden)) 1461 do isp=1,nspden 1462 do ii=1,qphase 1463 jj=(ii-1)*cplex*lmn2_size 1464 pawrhoij_gathered(jrhoij)%rhoijp(jj+1:jj+cplex*nselect,isp)= & 1465 & buf_dp_all(indx_dp:indx_dp+cplex*nselect-1) 1466 if (nselect<lmn2_size) & 1467 & pawrhoij_gathered(jrhoij)%rhoijp(jj+cplex*nselect+1:jj+cplex*lmn2_size,isp)=zero 1468 indx_dp=indx_dp+cplex*nselect 1469 end do 1470 end do 1471 end if 1472 if (lmnmix>0) then 1473 LIBPAW_ALLOCATE(pawrhoij_gathered(jrhoij)%kpawmix,(lmnmix)) 1474 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 if (ngrhoij>0) then 1478 LIBPAW_ALLOCATE(pawrhoij_gathered(jrhoij)%grhoij,(ngrhoij,qphase*cplex*lmn2_size,nspden)) 1479 do isp=1,nspden 1480 do ii=1,cplex*qphase*lmn2_size 1481 pawrhoij_gathered(jrhoij)%grhoij(1:ngrhoij,ii,isp)=buf_dp_all(indx_dp:indx_dp+ngrhoij-1) 1482 indx_dp=indx_dp+ngrhoij 1483 end do 1484 end do 1485 end if 1486 if (use_rhoijres>0) then 1487 LIBPAW_ALLOCATE(pawrhoij_gathered(jrhoij)%rhoijres,(qphase*cplex*lmn2_size,nspden)) 1488 do isp=1,nspden 1489 pawrhoij_gathered(jrhoij)%rhoijres(1:cplex*qphase*lmn2_size,isp)= & 1490 & buf_dp_all(indx_dp:indx_dp+cplex*qphase*lmn2_size-1) 1491 indx_dp=indx_dp+cplex*qphase*lmn2_size 1492 end do 1493 end if 1494 if (use_rhoij_>0) then 1495 LIBPAW_ALLOCATE(pawrhoij_gathered(jrhoij)%rhoij_,(qphase*cplex*lmn2_size,rhoij_size2)) 1496 do isp=1,rhoij_size2 1497 pawrhoij_gathered(jrhoij)%rhoij_(1:cplex*qphase*lmn2_size,isp)= & 1498 & buf_dp_all(indx_dp:indx_dp+cplex*qphase*lmn2_size-1) 1499 indx_dp=indx_dp+cplex*qphase*lmn2_size 1500 end do 1501 end if 1502 end do 1503 if ((indx_int/=1+buf_int_size_all).or.(indx_dp/=1+buf_dp_size_all)) then 1504 write(msg,*) 'Wrong buffer sizes: buf_int_size_all=',buf_int_size_all,' buf_dp_size_all=',buf_dp_size_all 1505 LIBPAW_BUG(msg) 1506 end if 1507 end if 1508 1509 !Free memory 1510 call free_my_atmtab(my_atmtab,my_atmtab_allocated) 1511 LIBPAW_DEALLOCATE(buf_int) 1512 LIBPAW_DEALLOCATE(buf_dp) 1513 LIBPAW_DEALLOCATE(buf_int_all) 1514 LIBPAW_DEALLOCATE(buf_dp_all) 1515 1516 end subroutine pawrhoij_gather
m_pawrhoij/pawrhoij_init_unpacked [ Functions ]
[ Top ] [ m_pawrhoij ] [ Functions ]
NAME
pawrhoij_init_unpacked
FUNCTION
Initialize field of rhoij datastructure for unpacked values (pawrhoij%rhoij_ array)
SIDE EFFECTS
rhoij(:) <pawrhoij_type)>= input/output datastructure * In output the rhoij_ array is allocated
SOURCE
2800 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 nrhoij=SIZE(rhoij);if (nrhoij==0) return 2813 2814 do iat=1,nrhoij 2815 2816 lmn2_size =rhoij(iat)%lmn2_size 2817 cplex_rhoij = rhoij(iat)%cplex_rhoij 2818 qphase = rhoij(iat)%qphase 2819 nsp2 = rhoij(iat)%nsppol;if (rhoij(iat)%nspden==4) nsp2=4 2820 2821 if (allocated(rhoij(iat)%rhoij_)) then 2822 LIBPAW_DEALLOCATE(rhoij(iat)%rhoij_) 2823 end if 2824 LIBPAW_ALLOCATE(rhoij(iat)%rhoij_,(cplex_rhoij*qphase*lmn2_size,nsp2)) 2825 rhoij(iat)%use_rhoij_=1 2826 rhoij(iat)%rhoij_=zero 2827 2828 end do 2829 2830 end subroutine pawrhoij_init_unpacked
m_pawrhoij/pawrhoij_inquire_dim [ Functions ]
[ Top ] [ m_pawrhoij ] [ Functions ]
NAME
pawrhoij_inquire_dim
FUNCTION
Compute the values f the dimensions (cplex_rhoij, qphase, nspden) for a pawrhoij datastructure. These ones depend on the parameters of the calculation
INPUTS
[cplex]= flag controlling the use of a exp(iqR) phase. 1=no phase, 2=phase [cpxocc]= 2 if rhoij is required to be imaginary [nspden]= number of spin-density components [qpt(3)]= q-vector, if any [spnorb]= flag: 1 if spin-orbit coupling is activated
OUTPUT
[cplex_rhoij] = value of cplex_rhoij associated to pawrhoij [qphase_rhoij]= value of qphase associated to pawrhoij [nspden_rhoij]= value of nspden associated to pawrhoij
SOURCE
3208 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 if (present(cplex_rhoij)) then 3225 cplex_rhoij=1 3226 if (present(cpxocc)) cplex_rhoij=max(cplex_rhoij,cpxocc) 3227 end if 3228 3229 !qphase_rhoij 3230 if (present(qphase_rhoij)) then 3231 qphase_rhoij=1 3232 if (present(cplex).and.present(qpt)) then 3233 msg='only one argument cplex or qpt should be passed!' 3234 LIBPAW_BUG(msg) 3235 end if 3236 if (present(cplex)) qphase_rhoij=merge(1,2,cplex==1) 3237 if (present(qpt)) then 3238 if (any(abs(qpt(:))>tol8)) qphase_rhoij=2 3239 end if 3240 end if 3241 3242 !nspden_rhoij 3243 if (present(nspden_rhoij)) then 3244 nspden_rhoij=1 3245 if (present(nspden)) nspden_rhoij=nspden 3246 if (present(spnorb)) nspden_rhoij=merge(nspden_rhoij,4,spnorb<=0) 3247 end if 3248 3249 end subroutine pawrhoij_inquire_dim
m_pawrhoij/pawrhoij_io [ Functions ]
[ Top ] [ m_pawrhoij ] [ Functions ]
NAME
pawrhoij_io
FUNCTION
IO method for pawrhoij datastructures.
INPUTS
unitfi=Unit number for IO file or netcdf file handler (already opened in the caller). nsppol_in=Number of independent spin polarizations. Only used for reading. nspinor_in=Number of spinorial components. Only used for reading. nspden_in=Number of spin-density components. only used for reading. nlmn_type(ntypat)= Number of (l,m,n) elements for the paw basis for each type of atom. Only used for reading. typat(natom) =Type of each atom. headform=Format of the abinit header (only used for reading as we need to know how to read the data. Writing is always done using the latest headform. rdwr_mode(len=*)=String defining the IO mode. Possible values (not case sensitive): "W"= For writing to unitfi "R"= For reading from unitfi "E"= For echoing. "D"= for debug [form(len=*)]= String defining the file format. Defaults to Fortran binary mode i.e., "unformatted" Other possible values are (case insensitive): "formatted"=For IO on a file open in formatted mode. "netcdf"=For IO on a netcdf file. [natinc]=Defines the increment in the loop over natom used for echoing the pawrhoij(natom) datastructures. If not specified, only the first and the last atom will be printed.
SIDE EFFECTS
pawrhoij(:)<type(pawrhoij_type)>= rhoij datastructure. if rdwr_mode="W", it will be written on unit unitfi using the file format given by form. if rdwr_mode="R", pawrhoij will be read and initialized from unit unitfi that has been opened with form=form. if rdwr_mode="E", the routines only echoes the content of the structure.
SOURCE
2280 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 integer,allocatable :: ibuffer(:),nsel44(:,:),nsel56(:) 2306 integer,pointer :: my_atmtab(:) 2307 real(dp), allocatable :: buffer(:) 2308 real(dp),pointer :: rhoij_tmp(:) 2309 2310 ! ************************************************************************* 2311 2312 my_natom=SIZE(pawrhoij);if (my_natom==0) return 2313 my_nspden=nspden_in 2314 natom=size(typat) 2315 paral_atom=(my_natom/=natom) 2316 if (present(mpi_atmtab)) then 2317 if (.not.associated(mpi_atmtab)) then 2318 msg='mpi_atmtab not associated (pawrhoij_io)' 2319 LIBPAW_BUG(msg) 2320 end if 2321 my_atmtab=>mpi_atmtab 2322 else if (my_natom/=natom) then 2323 msg='my_natom /=natom, mpi_atmtab should be in argument (pawrhoij_io)' 2324 LIBPAW_BUG(msg) 2325 end if 2326 2327 iomode = fort_binary 2328 if (PRESENT(form)) then 2329 select case (libpaw_to_upper(form)) 2330 case ("FORMATTED") 2331 iomode = fort_formatted 2332 case ("NETCDF") 2333 iomode = netcdf_io 2334 case default 2335 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 ncid = unitfi 2345 2346 select case (rdwr_mode(1:1)) 2347 2348 case ("R","r") ! Reading the Rhoij tab 2349 2350 if ((headform>=44).and.(headform<56)) then 2351 LIBPAW_ALLOCATE(nsel44,(nspden_in,natom)) 2352 if (iomode == fort_binary) then 2353 read(unitfi ) ((nsel44(ispden,iatom),ispden=1,nspden_in),iatom=1,natom) 2354 else if (iomode == fort_formatted) then 2355 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 LIBPAW_ERROR("header in 44-56 not compatible with Netcdf") 2359 #endif 2360 end if 2361 call pawrhoij_alloc(pawrhoij,1,nspden_in,nspinor_in,nsppol_in,typat,lmnsize=nlmn_type) 2362 do iatom=1,natom 2363 pawrhoij(iatom)%nrhoijsel=nsel44(1,iatom) 2364 end do 2365 bsize=sum(nsel44) 2366 LIBPAW_ALLOCATE(ibuffer,(bsize)) 2367 LIBPAW_ALLOCATE(buffer,(bsize)) 2368 if (iomode == fort_binary) then 2369 read(unitfi ) ibuffer(:),buffer(:) 2370 else if (iomode == fort_formatted) then 2371 read(unitfi,*) ibuffer(:),buffer(:) 2372 end if 2373 ii=0 2374 do iatom=1,natom 2375 nselect=nsel44(1,iatom) 2376 pawrhoij(iatom)%rhoijselect(:)=0 2377 pawrhoij(iatom)%rhoijselect(1:nselect)=ibuffer(ii+1:ii+nselect) 2378 do ispden=1,nspden_in 2379 pawrhoij(iatom)%rhoijp(1:nselect,ispden)=buffer(ii+1:ii+nselect) 2380 ii=ii+nselect 2381 end do 2382 end do 2383 LIBPAW_DEALLOCATE(ibuffer) 2384 LIBPAW_DEALLOCATE(buffer) 2385 LIBPAW_DEALLOCATE(nsel44) 2386 else if (headform>=56) then 2387 LIBPAW_ALLOCATE(nsel56,(natom)) 2388 my_cplex=1;my_nspden=1;my_qphase=1 2389 if (headform==56) then 2390 if (iomode == fort_binary) then 2391 read(unitfi ) (nsel56(iatom),iatom=1,natom),my_cplex 2392 else if (iomode == fort_formatted) then 2393 read(unitfi,*) (nsel56(iatom),iatom=1,natom),my_cplex 2394 #ifdef LIBPAW_HAVE_NETCDF 2395 else if (iomode == netcdf_io) then 2396 NCF_CHECK(nf90_inq_dimid(ncid, "pawrhoij_cplex", cplex_id)) 2397 NCF_CHECK(nf90_inquire_dimension(ncid, cplex_id, len=my_cplex)) 2398 2399 NCF_CHECK(nf90_inq_varid(ncid, "rhoijsel_atoms", nsel56_id)) 2400 NCF_CHECK(nf90_get_var(ncid, nsel56_id, nsel56)) 2401 #endif 2402 end if 2403 else 2404 if (iomode == fort_binary) then 2405 read(unitfi,err=10,end=10) (nsel56(iatom),iatom=1,natom),my_cplex,my_nspden,my_qphase 2406 10 continue 2407 else if (iomode == fort_formatted) then 2408 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 NCF_CHECK(nf90_inq_dimid(ncid, "pawrhoij_cplex", cplex_id)) 2413 NCF_CHECK(nf90_inquire_dimension(ncid, cplex_id, len=my_cplex)) 2414 NCF_CHECK(nf90_inq_dimid(ncid, "pawrhoij_nspden", nspden_id)) 2415 NCF_CHECK(nf90_inquire_dimension(ncid, nspden_id, len=my_nspden)) 2416 NCF_CHECK(nf90_inq_varid(ncid, "nrhoijsel_atoms", nsel56_id)) 2417 NCF_CHECK(nf90_get_var(ncid, nsel56_id, nsel56)) 2418 if (nf90_inq_dimid(ncid, "pawrhoij_qphase", qphase_id)==NF90_NOERR) then 2419 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 & lmnsize=nlmn_type,qphase=my_qphase) 2428 do iatom=1,natom 2429 pawrhoij(iatom)%nrhoijsel=nsel56(iatom) 2430 end do 2431 bsize=sum(nsel56) 2432 LIBPAW_ALLOCATE(ibuffer,(bsize)) 2433 LIBPAW_ALLOCATE(buffer,(bsize*my_nspden*my_cplex*my_qphase)) 2434 if (iomode == fort_binary) then 2435 read(unitfi ) ibuffer(:),buffer(:) 2436 else if (iomode == fort_formatted) then 2437 read(unitfi,*) ibuffer(:),buffer(:) 2438 #ifdef LIBPAW_HAVE_NETCDF 2439 else if (iomode == netcdf_io) then 2440 if (bsize > 0) then 2441 NCF_CHECK(nf90_inq_varid(ncid, "rhoijselect_atoms", ibuffer_id)) 2442 NCF_CHECK(nf90_get_var(ncid, ibuffer_id, ibuffer)) 2443 NCF_CHECK(nf90_inq_varid(ncid, "rhoijp_atoms", buffer_id)) 2444 NCF_CHECK(nf90_get_var(ncid, buffer_id, buffer)) 2445 end if 2446 #endif 2447 end if 2448 ii=0;jj=0 2449 do iatom=1,natom 2450 nselect=nsel56(iatom) 2451 pawrhoij(iatom)%rhoijselect(:)=0 2452 pawrhoij(iatom)%rhoijselect(1:nselect)=ibuffer(ii+1:ii+nselect) 2453 ii=ii+nselect 2454 do ispden=1,my_nspden 2455 pawrhoij(iatom)%rhoijp(1:my_cplex*nselect,ispden)= & 2456 buffer(jj+1:jj+my_cplex*nselect) 2457 jj=jj+my_cplex*nselect 2458 if (my_qphase==2) then 2459 itypat=typat(iatom) 2460 lmn_size=nlmn_type(itypat) 2461 lmn2_size=lmn_size*(lmn_size+1)/2 2462 iq0 = my_cplex*lmn2_size 2463 pawrhoij(iatom)%rhoijp(iq0+1:iq0+my_cplex*nselect,ispden)= & 2464 buffer(jj+1:jj+my_cplex*nselect) 2465 jj=jj+my_cplex*nselect 2466 end if 2467 end do 2468 end do 2469 LIBPAW_DEALLOCATE(ibuffer) 2470 LIBPAW_DEALLOCATE(buffer) 2471 LIBPAW_DEALLOCATE(nsel56) 2472 end if 2473 2474 case ("W","w") ! Writing the Rhoij tab (latest format is used) 2475 2476 LIBPAW_ALLOCATE(nsel56,(natom)) 2477 my_cplex =pawrhoij(1)%cplex_rhoij 2478 my_nspden=pawrhoij(1)%nspden 2479 my_qphase=pawrhoij(1)%qphase 2480 do iatom=1,natom 2481 nsel56(iatom)=pawrhoij(iatom)%nrhoijsel 2482 end do 2483 bsize=sum(nsel56) 2484 2485 if (iomode == fort_binary) then 2486 write(unitfi ) (nsel56(iatom),iatom=1,natom),my_cplex,my_nspden,my_qphase 2487 else if (iomode == fort_formatted) then 2488 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 ncerr = nf90_redef(ncid) 2492 2493 ! Define dimensions. 2494 ncerr = nf90_inq_dimid(ncid, "number_of_atoms", natom_id) 2495 if (ncerr /= nf90_noerr) then 2496 NCF_CHECK(nf90_def_dim(ncid, "number_of_atoms", natom, natom_id)) 2497 end if 2498 ncerr = nf90_inq_varid(ncid, "nrhoijsel_atoms", nsel56_id) 2499 if (ncerr /= nf90_noerr) then 2500 NCF_CHECK(nf90_def_var(ncid, "nrhoijsel_atoms", NF90_INT, natom_id, nsel56_id)) 2501 end if 2502 ncerr = nf90_inq_dimid(ncid, "pawrhoij_cplex", cplex_id) 2503 if (ncerr /= nf90_noerr) then 2504 NCF_CHECK(nf90_def_dim(ncid, "pawrhoij_cplex", my_cplex, cplex_id)) 2505 end if 2506 ncerr = nf90_inq_dimid(ncid, "pawrhoij_nspden", nspden_id) 2507 if (ncerr /= nf90_noerr) then 2508 NCF_CHECK(nf90_def_dim(ncid, "pawrhoij_nspden", my_nspden, nspden_id)) 2509 end if 2510 ncerr = nf90_inq_dimid(ncid, "pawrhoij_qphase", qphase_id) 2511 if (ncerr /= nf90_noerr) then 2512 NCF_CHECK(nf90_def_dim(ncid, "pawrhoij_qphase", my_qphase, qphase_id)) 2513 end if 2514 if (bsize > 0) then 2515 ncerr = nf90_inq_dimid(ncid, "rhoijselect_atoms_dim", bsize_id) 2516 if (ncerr /= nf90_noerr) then 2517 NCF_CHECK(nf90_def_dim(ncid, "rhoijselect_atoms_dim", bsize, bsize_id)) 2518 end if 2519 ncerr = nf90_inq_dimid(ncid, "rhoijp_atoms_dim", bufsize_id) 2520 if (ncerr /= nf90_noerr) then 2521 NCF_CHECK(nf90_def_dim(ncid, "rhoijp_atoms_dim", bsize*my_nspden*my_qphase*my_cplex, bufsize_id)) 2522 end if 2523 ncerr = nf90_inq_varid(ncid, "rhoijselect_atoms", ibuffer_id) 2524 if (ncerr /= nf90_noerr) then 2525 NCF_CHECK(nf90_def_var(ncid, "rhoijselect_atoms", NF90_INT, bsize_id, ibuffer_id)) 2526 end if 2527 ncerr = nf90_inq_varid(ncid, "rhoijp_atoms", buffer_id) 2528 if (ncerr /= nf90_noerr) then 2529 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 LIBPAW_COMMENT("All rhoij entries are zero. No netcdf entry produced") 2534 end if 2535 2536 ! Write nsel56 2537 NCF_CHECK(nf90_enddef(ncid)) 2538 NCF_CHECK(nf90_put_var(ncid, nsel56_id, nsel56)) 2539 #endif 2540 end if 2541 2542 LIBPAW_ALLOCATE(ibuffer,(bsize)) 2543 LIBPAW_ALLOCATE(buffer,(bsize*my_nspden*my_cplex*my_qphase)) 2544 ii=0;jj=0 2545 do iatom=1,natom 2546 nselect=nsel56(iatom) 2547 ibuffer(ii+1:ii+nselect)=pawrhoij(iatom)%rhoijselect(1:nselect) 2548 ii=ii+nselect 2549 do ispden=1,my_nspden 2550 buffer(jj+1:jj+my_cplex*nselect)=& 2551 pawrhoij(iatom)%rhoijp(1:my_cplex*nselect,ispden) 2552 jj=jj+my_cplex*nselect 2553 if (my_qphase==2) then 2554 iq0 = my_cplex*pawrhoij(iatom)%lmn2_size 2555 buffer(jj+1:jj+my_cplex*nselect)=& 2556 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 if (iomode == fort_binary) then 2562 write(unitfi ) ibuffer(:),buffer(:) 2563 else if (iomode == fort_formatted) then 2564 write(unitfi,*) ibuffer(:),buffer(:) 2565 #ifdef LIBPAW_HAVE_NETCDF 2566 else if (iomode == netcdf_io) then 2567 if (bsize > 0) then 2568 NCF_CHECK(nf90_put_var(ncid, ibuffer_id, ibuffer)) 2569 NCF_CHECK(nf90_put_var(ncid, buffer_id, buffer)) 2570 end if 2571 #endif 2572 end if 2573 LIBPAW_DEALLOCATE(ibuffer) 2574 LIBPAW_DEALLOCATE(buffer) 2575 LIBPAW_DEALLOCATE(nsel56) 2576 2577 case ("E","e") ! Echoing the Rhoij tab 2578 2579 my_natinc=1; if(natom>1) my_natinc=natom-1 2580 my_qphase=pawrhoij(1)%qphase 2581 nselect=maxval(pawrhoij(:)%nrhoijsel) 2582 if (PRESENT(natinc)) my_natinc = natinc ! user-defined increment. 2583 LIBPAW_ALLOCATE(ibuffer,(0)) 2584 nselect=maxval(pawrhoij(:)%nrhoijsel) 2585 if (my_qphase==2) then 2586 LIBPAW_POINTER_ALLOCATE(rhoij_tmp,(2*nselect)) 2587 end if 2588 do iatom=1,my_natom,my_natinc 2589 iatom_tot=iatom;if(paral_atom)iatom_tot=my_atmtab(iatom) 2590 my_cplex =pawrhoij(iatom)%cplex_rhoij 2591 my_nspden=pawrhoij(iatom)%nspden 2592 my_qphase=pawrhoij(iatom)%qphase 2593 nselect=pawrhoij(iatom)%nrhoijsel 2594 do ispden=1,pawrhoij(iatom)%nspden 2595 if (my_qphase==1) then 2596 my_cplex_eff=my_cplex 2597 rhoij_tmp => pawrhoij(iatom)%rhoijp(1:my_cplex*nselect,ispden) 2598 else 2599 my_cplex_eff=2 2600 rhoij_tmp=zero 2601 jj=my_cplex*pawrhoij(iatom)%lmn2_size 2602 if (my_cplex==1) then 2603 do ii=1,nselect 2604 rhoij_tmp(2*ii-1)=pawrhoij(iatom)%rhoijp(ii,ispden) 2605 rhoij_tmp(2*ii )=pawrhoij(iatom)%rhoijp(jj+ii,ispden) 2606 end do 2607 else 2608 do ii=1,nselect 2609 rhoij_tmp(2*ii-1)=pawrhoij(iatom)%rhoijp(2*ii-1,ispden) & 2610 & -pawrhoij(iatom)%rhoijp(jj+2*ii ,ispden) 2611 rhoij_tmp(2*ii )=pawrhoij(iatom)%rhoijp(2*ii ,ispden) & 2612 & +pawrhoij(iatom)%rhoijp(jj+2*ii-1,ispden) 2613 end do 2614 end if 2615 end if 2616 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 & opt_sym=2,mode_paral='PERS') 2621 end do ! end nspden do 2622 end do ! end iatom do 2623 LIBPAW_DEALLOCATE(ibuffer) 2624 if (my_qphase==2) then 2625 LIBPAW_POINTER_DEALLOCATE(rhoij_tmp) 2626 end if 2627 2628 case ("D","d") ! Debug 2629 2630 write(unitfi,'(a,i4)' ) 'size pawmkrhoij , natom = ' , my_natom 2631 my_natinc=1; if(natom>1) my_natinc=natom-1 2632 if (PRESENT(natinc)) my_natinc = natinc ! user-defined increment. 2633 LIBPAW_ALLOCATE(ibuffer,(0)) 2634 do iatom=1,my_natom,my_natinc 2635 iatom_tot=iatom;if(paral_atom) iatom_tot=my_atmtab(iatom) 2636 if (iatom_tot/=1) cycle 2637 write(unitfi,'(a,i4,a)' ) ' ******* rhoij (Atom # ',iatom_tot,' ********)' 2638 write(unitfi,'(a,i4,a,i4)' ) 'cplex_rhoij=',pawrhoij(iatom)%cplex_rhoij, ' nselect=', pawrhoij(iatom)%nrhoijsel 2639 write(unitfi,'(a,i4,a,i4)' ) 'nspden=',pawrhoij(iatom)%nspden, ' lmn2size=',pawrhoij(iatom)%lmn2_size 2640 write(unitfi,'(a,i4,a,i4)' ) 'lmnmix=',pawrhoij(iatom)%lmnmix_sz, ' ngrhoij=',pawrhoij(iatom)%ngrhoij 2641 write(unitfi,'(a,i4,a,i4)' ) 'use_rhoijres=',pawrhoij(iatom)%use_rhoijres, & 2642 & 'use_rhoij_=',pawrhoij(iatom)%use_rhoij_ 2643 write(unitfi,'(a,i4,a,i4)' ) 'itypat=',pawrhoij(iatom)%itypat, ' lmn_size=',pawrhoij(iatom)%lmn_size 2644 write(unitfi,'(a,i4,a,i4)' ) 'nsppol=',pawrhoij(iatom)%nsppol, ' nspinor=',pawrhoij(iatom)%nspinor 2645 write(unitfi,'(a,i4)' ) 'qphase=',pawrhoij(iatom)%qphase 2646 cplex=pawrhoij(iatom)%cplex_rhoij 2647 qphase=pawrhoij(iatom)%qphase 2648 lmn2_size=pawrhoij(iatom)%lmn2_size 2649 do i2=1,pawrhoij(iatom)%nrhoijsel 2650 write(unitfi,'(a,i4,a,i4,a,i4,a,f9.5)') 'rhoijselect(,',i2,')=',& 2651 & pawrhoij(iatom)%rhoijselect(i2) 2652 end do 2653 if (pawrhoij(iatom)%ngrhoij>0) then 2654 ngrhoijmx=2; if (pawrhoij(iatom)%ngrhoij<ngrhoijmx) ngrhoijmx=pawrhoij(iatom)%ngrhoij 2655 do ispden=1,pawrhoij(iatom)%nspden 2656 do i1=ngrhoijmx,ngrhoijmx 2657 do ii=1,qphase 2658 do i2=(ii-1)*cplex*lmn2_size+cplex*lmn2_size,(ii-1)*cplex*lmn2_size+cplex*lmn2_size 2659 write(unitfi,'(a,i4,a,i4,a,i4,a,f9.5)') ' grhoij(,',i1,',',i2,',',ispden,')=',& 2660 & pawrhoij(iatom)%grhoij(i1,i2,ispden) 2661 end do 2662 end do 2663 end do 2664 end do 2665 call libpaw_flush(unitfi) 2666 end if 2667 if (pawrhoij(iatom)%use_rhoijres>0) then 2668 do ispden=1,pawrhoij(iatom)%nspden 2669 do ii=1,qphase 2670 do i2=(ii-1)*cplex*lmn2_size+cplex*lmn2_size,(ii-1)*cplex*lmn2_size+cplex*lmn2_size 2671 write(unitfi,'(a,i4,a,i4,a,f9.5)') ' rhoijres(,',i2,',ispden=',ispden,')=',& 2672 & pawrhoij(iatom)%rhoijres(i2,ispden) 2673 end do 2674 end do 2675 end do 2676 call libpaw_flush(unitfi) 2677 end if 2678 if (pawrhoij(iatom)%nrhoijsel>0) then 2679 do ispden=1,pawrhoij(iatom)%nspden 2680 do ii=1,qphase 2681 do i2=(ii-1)*cplex*lmn2_size+cplex*pawrhoij(iatom)%nrhoijsel, & 2682 & (ii-1)*cplex*lmn2_size+cplex*pawrhoij(iatom)%nrhoijsel 2683 write(unitfi,'(a,i4,a,i4,a,f9.5)') ' rhoijp!(nrhoijselec=,',i2,',ispden=',ispden,')=',& 2684 & pawrhoij(iatom)%rhoijp(i2,ispden) 2685 end do 2686 end do 2687 end do 2688 call libpaw_flush(unitfi) 2689 end if 2690 if (pawrhoij(iatom)%use_rhoij_>0) then 2691 size_rhoij2=size(pawrhoij(iatom)%rhoij_,2) 2692 do ispden=1,size_rhoij2 2693 do ii=1,qphase 2694 do i2=(ii-1)*cplex*lmn2_size+cplex*lmn2_size,(ii-1)*cplex*lmn2_size+cplex*lmn2_size 2695 write(unitfi,'(a,i4,a,i4,a,f9.5)') ' rhoij_(,',i2,',ispden=',ispden,')=',& 2696 & pawrhoij(iatom)%rhoij_(i2,ispden) 2697 end do 2698 end do 2699 end do 2700 end if 2701 call libpaw_flush(unitfi) 2702 if (pawrhoij(iatom)%lmnmix_sz>0) then 2703 write(unitfi,'(a)') 'kpawmix ' 2704 write(unitfi,'(i4,i4,i4,i4)') pawrhoij(iatom)%kpawmix(:) 2705 end if 2706 call libpaw_flush(unitfi) 2707 end do 2708 2709 case default 2710 msg='Wrong rdwr_mode'//TRIM(rdwr_mode) 2711 LIBPAW_ERROR(msg) 2712 2713 end select 2714 2715 end subroutine pawrhoij_io
m_pawrhoij/pawrhoij_isendreceive_fillbuffer [ Functions ]
[ Top ] [ m_pawrhoij ] [ Functions ]
NAME
pawrhoij_isendreceive_fillbuffer
FUNCTION
Extract from pawrhoij and from the global index of atoms the buffers to send in a sending operation This function has to be coupled with a call to pawrhoij_isendreceive_getbuffer
INPUTS
atm_indx_send(1:total number of atoms)= array for send operation,
Given an index of atom in global numbering, give its index
in the table of atoms treated by current processor
or -1 if the atoms is not treated by current processor
nrhoij_send= number of sent atoms
pawrhoij= data structure from which are extract buffer int and buffer dp
OUTPUT
buf_int= buffer of integers to be sent buf_int_size= size of buffer of integers buf_dp= buffer of double precision numbers to be sent buf_dp_size= size of buffer of double precision numbers
SOURCE
4487 subroutine pawrhoij_isendreceive_fillbuffer(pawrhoij,atmtab_send, atm_indx_send,nrhoij_send,& 4488 & buf_int,buf_int_size,buf_dp,buf_dp_size) 4489 4490 !Arguments ------------------------------------ 4491 !scalars 4492 integer,intent(out) :: buf_dp_size,buf_int_size 4493 integer,intent(in) :: nrhoij_send 4494 !arrays 4495 integer,intent(in) :: atmtab_send(:),atm_indx_send(:) 4496 integer,intent(out),allocatable :: buf_int(:) 4497 real(dp),intent(out),allocatable :: buf_dp(:) 4498 type(pawrhoij_type),target,intent(in) :: pawrhoij(:) 4499 4500 !Local variables------------------------------- 4501 !scalars 4502 integer :: cplex,ii,indx_int,indx_dp, iatom_tot,irhoij,irhoij_send,isp,jj,lmn2_size,lmnmix 4503 integer :: ngrhoij,nselect,nspden,qphase,rhoij_size2 4504 integer :: use_rhoijp,use_rhoijres,use_rhoij_ 4505 character(len=500) :: msg 4506 type(pawrhoij_type),pointer :: pawrhoij1 4507 !arrays 4508 4509 ! ********************************************************************* 4510 4511 !Compute sizes of buffers 4512 buf_int_size=0;buf_dp_size=0 4513 nselect=0;lmnmix=0;ngrhoij=0;rhoij_size2=0 4514 use_rhoijp=0;use_rhoijres=0;use_rhoij_=0 4515 do irhoij_send=1,nrhoij_send 4516 iatom_tot=atmtab_send(irhoij_send) 4517 irhoij=atm_indx_send(iatom_tot) 4518 if (irhoij == -1) then 4519 msg="Error in pawrhoij_isendreceive_fillbuffer atom not found" 4520 LIBPAW_BUG(msg) 4521 end if 4522 pawrhoij1=>pawrhoij(irhoij) 4523 cplex =pawrhoij1%cplex_rhoij 4524 qphase =pawrhoij1%qphase 4525 lmn2_size=pawrhoij1%lmn2_size 4526 nspden =pawrhoij1%nspden 4527 lmnmix=pawrhoij1%lmnmix_sz 4528 ngrhoij=pawrhoij1%ngrhoij 4529 use_rhoijp=pawrhoij1%use_rhoijp 4530 use_rhoijres=pawrhoij1%use_rhoijres 4531 use_rhoij_=pawrhoij1%use_rhoij_ 4532 buf_int_size=buf_int_size+16 4533 if (use_rhoijp>0) then 4534 nselect=pawrhoij1%nrhoijsel 4535 buf_int_size=buf_int_size+nselect 4536 buf_dp_size=buf_dp_size + cplex*qphase*nselect*nspden 4537 end if 4538 if (lmnmix>0) buf_int_size=buf_int_size+lmnmix 4539 if (ngrhoij>0) buf_dp_size=buf_dp_size + cplex*qphase*lmn2_size*nspden*ngrhoij 4540 if (use_rhoijres>0) buf_dp_size=buf_dp_size + cplex*qphase*lmn2_size*nspden 4541 if (use_rhoij_>0) then 4542 rhoij_size2=size(pawrhoij1%rhoij_,dim=2) 4543 buf_dp_size=buf_dp_size + cplex*qphase*lmn2_size*rhoij_size2 4544 end if 4545 end do 4546 4547 !Fill input buffers 4548 LIBPAW_ALLOCATE(buf_int,(buf_int_size)) 4549 LIBPAW_ALLOCATE(buf_dp,(buf_dp_size)) 4550 indx_int=1;indx_dp =1 4551 lmnmix=0;ngrhoij=0;nselect=0;rhoij_size2=0 4552 use_rhoijp=0;use_rhoijres=0;use_rhoij_=0 4553 do irhoij_send=1,nrhoij_send 4554 iatom_tot=atmtab_send(irhoij_send) 4555 irhoij=atm_indx_send(iatom_tot) 4556 pawrhoij1=>pawrhoij(irhoij) 4557 cplex =pawrhoij1%cplex_rhoij 4558 qphase =pawrhoij1%qphase 4559 lmn2_size=pawrhoij1%lmn2_size 4560 nspden =pawrhoij1%nspden 4561 lmnmix=pawrhoij1%lmnmix_sz 4562 ngrhoij=pawrhoij1%ngrhoij 4563 use_rhoijp=pawrhoij1%use_rhoijp 4564 nselect=pawrhoij1%nrhoijsel 4565 use_rhoijres=pawrhoij1%use_rhoijres 4566 use_rhoij_ =pawrhoij1%use_rhoij_ 4567 if (use_rhoij_ > 0) then 4568 rhoij_size2 =size(pawrhoij1%rhoij_,dim=2) 4569 end if 4570 buf_int(indx_int)=atmtab_send(irhoij_send) ;indx_int=indx_int+1 4571 buf_int(indx_int)=cplex ;indx_int=indx_int+1 4572 buf_int(indx_int)=qphase ;indx_int=indx_int+1 4573 buf_int(indx_int)=lmn2_size ;indx_int=indx_int+1 4574 buf_int(indx_int)=nspden ;indx_int=indx_int+1 4575 buf_int(indx_int)=nselect ;indx_int=indx_int+1 4576 buf_int(indx_int)=lmnmix ;indx_int=indx_int+1 4577 buf_int(indx_int)=ngrhoij ;indx_int=indx_int+1 4578 buf_int(indx_int)=use_rhoijp ;indx_int=indx_int+1 4579 buf_int(indx_int)=use_rhoijres ;indx_int=indx_int+1 4580 buf_int(indx_int)=use_rhoij_ ;indx_int=indx_int+1 4581 buf_int(indx_int)=rhoij_size2 ;indx_int=indx_int+1 4582 buf_int(indx_int)=pawrhoij1%itypat ;indx_int=indx_int+1 4583 buf_int(indx_int)=pawrhoij1%lmn_size ;indx_int=indx_int+1 4584 buf_int(indx_int)=pawrhoij1%nsppol ;indx_int=indx_int+1 4585 buf_int(indx_int)=pawrhoij1%nspinor ;indx_int=indx_int+1 4586 if (use_rhoijp>0) then 4587 buf_int(indx_int:indx_int+nselect-1)=pawrhoij1%rhoijselect(1:nselect) 4588 indx_int=indx_int+nselect 4589 do isp=1,nspden 4590 do ii=1,qphase 4591 jj=(ii-1)*cplex*lmn2_size 4592 buf_dp(indx_dp:indx_dp+cplex*nselect-1)=pawrhoij1%rhoijp(jj+1:jj+cplex*nselect,isp) 4593 indx_dp=indx_dp+cplex*nselect 4594 end do 4595 end do 4596 end if 4597 if (lmnmix>0) then 4598 buf_int(indx_int:indx_int+lmnmix-1)=pawrhoij1%kpawmix(1:lmnmix) 4599 indx_int=indx_int+lmnmix 4600 end if 4601 if (ngrhoij>0) then 4602 do isp=1,nspden 4603 do ii=1,cplex*qphase*lmn2_size 4604 buf_dp(indx_dp:indx_dp+ngrhoij-1)=pawrhoij1%grhoij(1:ngrhoij,ii,isp) 4605 indx_dp=indx_dp+ngrhoij 4606 end do 4607 end do 4608 end if 4609 if (use_rhoijres>0) then 4610 do isp=1,nspden 4611 buf_dp(indx_dp:indx_dp+cplex*qphase*lmn2_size-1)=pawrhoij1%rhoijres(1:cplex*qphase*lmn2_size,isp) 4612 indx_dp=indx_dp+cplex*qphase*lmn2_size 4613 end do 4614 end if 4615 if (use_rhoij_>0) then 4616 do isp=1,rhoij_size2 4617 buf_dp(indx_dp:indx_dp+cplex*qphase*lmn2_size-1)=pawrhoij1%rhoij_(1:cplex*qphase*lmn2_size,isp) 4618 indx_dp=indx_dp+cplex*qphase*lmn2_size 4619 end do 4620 end if 4621 end do !irhoij_send 4622 4623 !Check 4624 if ((indx_int-1/=buf_int_size).or.(indx_dp-1/=buf_dp_size)) then 4625 write(msg,'(a,i10,a,i10)') 'Wrong buffer sizes: buf_int_size=',buf_int_size,' buf_dp_size=',buf_dp_size 4626 LIBPAW_BUG(msg) 4627 end if 4628 4629 end subroutine pawrhoij_isendreceive_fillbuffer
m_pawrhoij/pawrhoij_isendreceive_getbuffer [ Functions ]
[ Top ] [ m_pawrhoij ] [ Functions ]
NAME
pawrhoij_isendreceive_getbuffer
FUNCTION
Fill a pawrhoij structure with the buffers received in a receive operation This buffer should have been first extracted by a call to pawrhoij_isendreceive_fillbuffer
INPUTS
atm_indx_recv(1:total number of atoms)= array for receive operation
Given an index of atom in global numbering, give its index
in the table of atoms treated by current processor
or -1 if the atoms is not treated by current processor
buf_int= buffer of receive integers
buf_dp= buffer of receive double precision numbers
nrhoij_send= number of sent atoms
OUTPUT
pawrhoij= output datastructure filled with buffers receive in a receive operation
SOURCE
4348 subroutine pawrhoij_isendreceive_getbuffer(pawrhoij,nrhoij_send,atm_indx_recv,buf_int,buf_dp) 4349 4350 !Arguments ------------------------------------ 4351 !scalars 4352 integer,intent(in) :: nrhoij_send 4353 !arrays 4354 integer,intent(in) ::atm_indx_recv(:),buf_int(:) 4355 real(dp),intent(in) :: buf_dp(:) 4356 type(pawrhoij_type),target,intent(inout) :: pawrhoij(:) 4357 4358 !Local variables------------------------------- 4359 !scalars 4360 integer :: buf_dp_size,buf_int_size,cplex,ii,indx_int,indx_dp,iatom_tot,irhoij_send 4361 integer :: isp,jj,jrhoij,lmn2_size,lmnmix,ngrhoij,nselect,nspden,qphase,rhoij_size2 4362 integer :: use_rhoijp,use_rhoijres,use_rhoij_ 4363 character(len=500) :: msg 4364 type(pawrhoij_type),pointer :: pawrhoij1 4365 !arrays 4366 4367 ! ********************************************************************* 4368 4369 buf_int_size=size(buf_int) 4370 buf_dp_size=size(buf_dp) 4371 indx_int=1;indx_dp=1 4372 4373 do irhoij_send=1,nrhoij_send 4374 4375 iatom_tot=buf_int(indx_int) ; indx_int=indx_int+1 4376 jrhoij= atm_indx_recv(iatom_tot) 4377 if (jrhoij==-1) then 4378 msg="Error in pawrhoij_isendreceive_getbuffer atom not found" 4379 LIBPAW_BUG(msg) 4380 end if 4381 pawrhoij1=>pawrhoij(jrhoij) 4382 4383 cplex =buf_int(indx_int) ;indx_int=indx_int+1 4384 qphase =buf_int(indx_int) ;indx_int=indx_int+1 4385 lmn2_size =buf_int(indx_int) ;indx_int=indx_int+1 4386 nspden =buf_int(indx_int) ;indx_int=indx_int+1 4387 nselect =buf_int(indx_int) ;indx_int=indx_int+1 4388 lmnmix =buf_int(indx_int) ;indx_int=indx_int+1 4389 ngrhoij =buf_int(indx_int) ;indx_int=indx_int+1 4390 use_rhoijp =buf_int(indx_int) ;indx_int=indx_int+1 4391 use_rhoijres=buf_int(indx_int) ;indx_int=indx_int+1 4392 use_rhoij_ =buf_int(indx_int) ;indx_int=indx_int+1 4393 rhoij_size2 =buf_int(indx_int) ;indx_int=indx_int+1 4394 pawrhoij1%itypat=buf_int(indx_int) ;indx_int=indx_int+1 4395 pawrhoij1%lmn_size=buf_int(indx_int) ;indx_int=indx_int+1 4396 pawrhoij1%nsppol=buf_int(indx_int) ;indx_int=indx_int+1 4397 pawrhoij1%nspinor=buf_int(indx_int) ;indx_int=indx_int+1 4398 pawrhoij1%cplex_rhoij=cplex 4399 pawrhoij1%qphase=qphase 4400 pawrhoij1%lmn2_size=lmn2_size 4401 pawrhoij1%nspden=nspden 4402 pawrhoij1%nrhoijsel=nselect 4403 pawrhoij1%lmnmix_sz=lmnmix 4404 pawrhoij1%ngrhoij=ngrhoij 4405 pawrhoij1%use_rhoijp=use_rhoijp 4406 pawrhoij1%use_rhoijres=use_rhoijres 4407 pawrhoij1%use_rhoij_=use_rhoij_ 4408 if (use_rhoijp>0) then 4409 LIBPAW_ALLOCATE(pawrhoij1%rhoijselect,(lmn2_size)) 4410 pawrhoij1%rhoijselect(1:nselect)=buf_int(indx_int:indx_int+nselect-1) 4411 if (nselect < lmn2_size )pawrhoij1%rhoijselect(nselect+1:lmn2_size)=zero 4412 indx_int=indx_int+nselect 4413 LIBPAW_ALLOCATE(pawrhoij1%rhoijp,(cplex*qphase*lmn2_size,nspden)) 4414 do isp=1,nspden 4415 do ii=1,qphase 4416 jj=(ii-1)*cplex*lmn2_size 4417 pawrhoij1%rhoijp(jj+1:jj+cplex*nselect,isp)=buf_dp(indx_dp:indx_dp+cplex*nselect-1) 4418 if (nselect<lmn2_size)pawrhoij1%rhoijp(jj+cplex*nselect+1:jj+cplex*lmn2_size,isp)=zero 4419 indx_dp=indx_dp+cplex*nselect 4420 end do 4421 end do 4422 end if 4423 if (lmnmix>0) then 4424 LIBPAW_ALLOCATE(pawrhoij1%kpawmix,(lmnmix)) 4425 pawrhoij1%kpawmix(1:lmnmix)=buf_int(indx_int:indx_int+lmnmix-1) 4426 indx_int=indx_int+lmnmix 4427 end if 4428 if (ngrhoij>0) then 4429 LIBPAW_ALLOCATE(pawrhoij1%grhoij,(ngrhoij,cplex*qphase*lmn2_size,nspden)) 4430 do isp=1,nspden 4431 do ii=1,cplex*qphase*lmn2_size 4432 pawrhoij1%grhoij(1:ngrhoij,ii,isp)=buf_dp(indx_dp:indx_dp+ngrhoij-1) 4433 indx_dp=indx_dp+ngrhoij 4434 end do 4435 end do 4436 end if 4437 if (use_rhoijres>0) then 4438 LIBPAW_ALLOCATE(pawrhoij1%rhoijres,(cplex*qphase*lmn2_size,nspden)) 4439 do isp=1,nspden 4440 pawrhoij1%rhoijres(1:cplex*qphase*lmn2_size,isp)=buf_dp(indx_dp:indx_dp+cplex*qphase*lmn2_size-1) 4441 indx_dp=indx_dp+cplex*qphase*lmn2_size 4442 end do 4443 end if 4444 if (use_rhoij_>0) then 4445 LIBPAW_ALLOCATE(pawrhoij1%rhoij_,(cplex*qphase*lmn2_size,rhoij_size2)) 4446 do isp=1,rhoij_size2 4447 pawrhoij1%rhoij_(1:cplex*qphase*lmn2_size,isp)=buf_dp(indx_dp:indx_dp+cplex*qphase*lmn2_size-1) 4448 indx_dp=indx_dp+cplex*qphase*lmn2_size 4449 end do 4450 end if 4451 end do !irhoij_send 4452 if ((indx_int/=1+buf_int_size).or.(indx_dp/=1+buf_dp_size)) then 4453 write(msg,'(a,i10,a,i10)') 'Wrong buffer sizes: buf_int_size=',buf_int_size,' buf_dp_size=',buf_dp_size 4454 LIBPAW_BUG(msg) 4455 end if 4456 4457 end subroutine pawrhoij_isendreceive_getbuffer
m_pawrhoij/pawrhoij_mpisum_unpacked_1D [ Functions ]
[ Top ] [ m_pawrhoij ] [ Functions ]
NAME
pawrhoij_mpisum_unpacked_1D
FUNCTION
Build the MPI sum of the unsymmetrized PAW rhoij_ (augmentation occupancies)
Remember:for each atom, rho_ij=Sum_{n,k} {occ(n,k)*<Cnk|p_i><p_j|Cnk>}
Target: 1D array of pawrhoij datastructures
INPUTS
comm1=MPI communicator. Data will be MPI summed inside comm1
[comm2]=second MPI communicator. If present, rhoij_ will be
MPI summed inside comm2 after the collective sum in comm1.
SIDE EFFECTS
pawrhoij(:) <type(pawrhoij_type)>= paw rhoij occupancies and related data Input: the data calculateed by this processor. Output: the final MPI sum over comm1 and comm2.
SOURCE
2897 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 integer,allocatable :: dimlmn(:) 2913 real(dp),allocatable :: buffer(:) 2914 2915 !************************************************************************ 2916 2917 natom=SIZE(pawrhoij);if (natom==0) return 2918 2919 nproc1 = xmpi_comm_size(comm1) 2920 nproc2=1; if (PRESENT(comm2)) nproc2 = xmpi_comm_size(comm2) 2921 if (nproc1==1.and.nproc2==1) RETURN 2922 2923 !Fill the MPI buffer from the local rhoij_ 2924 LIBPAW_ALLOCATE(dimlmn,(natom)) 2925 dimlmn(1:natom)=pawrhoij(1:natom)%cplex_rhoij & 2926 & *pawrhoij(1:natom)%qphase & 2927 & *pawrhoij(1:natom)%lmn2_size 2928 nsp2=pawrhoij(1)%nsppol; if (pawrhoij(1)%nspden==4) nsp2=4 2929 bufdim=sum(dimlmn)*nsp2 2930 LIBPAW_ALLOCATE(buffer,(bufdim)) 2931 jdim=0 2932 do iatom=1,natom 2933 do isppol=1,nsp2 2934 buffer(jdim+1:jdim+dimlmn(iatom))=pawrhoij(iatom)%rhoij_(:,isppol) 2935 jdim=jdim+dimlmn(iatom) 2936 end do 2937 end do 2938 2939 !Build sum of pawrhoij%rhoij_ 2940 call xmpi_sum(buffer,comm1,ierr) ! Sum over the first communicator. 2941 if (PRESENT(comm2)) then 2942 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 do iatom=1,natom 2948 do isppol=1,nsp2 2949 pawrhoij(iatom)%rhoij_(:,isppol)=buffer(jdim+1:jdim+dimlmn(iatom)) 2950 jdim=jdim+dimlmn(iatom) 2951 end do 2952 end do 2953 2954 LIBPAW_DEALLOCATE(buffer) 2955 LIBPAW_DEALLOCATE(dimlmn) 2956 2957 end subroutine pawrhoij_mpisum_unpacked_1D
m_pawrhoij/pawrhoij_mpisum_unpacked_2D [ Functions ]
[ Top ] [ m_pawrhoij ] [ Functions ]
NAME
pawrhoij_mpisum_unpacked_2D
FUNCTION
Build the MPI sum of the unsymmetrized PAW rhoij_ (augmentation occupancies)
Remember:for each atom, rho_ij=Sum_{n,k} {occ(n,k)*<Cnk|p_i><p_j|Cnk>}
Target: 2D array of pawrhoij datastructures
INPUTS
comm1=MPI communicator. Data will be MPI summed inside comm1
[comm2]=second MPI communicator. If present, rhoij_ will be
MPI summed inside comm2 after the collective sum in comm1.
SIDE EFFECTS
pawrhoij(:,:) <type(pawrhoij_type)>= paw rhoij occupancies and related data Input: the data calculateed by this processor. Output: the final MPI sum over comm1 and comm2.
SOURCE
2983 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 integer,allocatable :: dimlmn(:,:) 2999 real(dp),allocatable :: buffer(:) 3000 3001 !************************************************************************ 3002 3003 natom =SIZE(pawrhoij,1);if (natom ==0) return 3004 nrhoij=SIZE(pawrhoij,2);if (nrhoij==0) return 3005 3006 nproc1 = xmpi_comm_size(comm1) 3007 nproc2=1; if (PRESENT(comm2)) nproc2 = xmpi_comm_size(comm2) 3008 if (nproc1==1.and.nproc2==1) RETURN 3009 3010 !Fill the MPI buffer from the local rhoij_ 3011 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 & *pawrhoij(1:natom,1:nrhoij)%lmn2_size 3015 nsp2=pawrhoij(1,1)%nsppol; if (pawrhoij(1,1)%nspden==4) nsp2=4 3016 bufdim=sum(dimlmn)*nsp2 3017 LIBPAW_ALLOCATE(buffer,(bufdim)) 3018 jdim=0 3019 do irhoij=1,nrhoij 3020 do iatom=1,natom 3021 do isppol=1,nsp2 3022 buffer(jdim+1:jdim+dimlmn(iatom,irhoij))=pawrhoij(iatom,irhoij)%rhoij_(:,isppol) 3023 jdim=jdim+dimlmn(iatom,irhoij) 3024 end do 3025 end do 3026 end do 3027 3028 !Build sum of pawrhoij%rhoij_ 3029 call xmpi_sum(buffer,comm1,ierr) ! Sum over the first communicator. 3030 if (PRESENT(comm2)) then 3031 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 do irhoij=1,nrhoij 3037 do iatom=1,natom 3038 do isppol=1,nsp2 3039 pawrhoij(iatom,irhoij)%rhoij_(:,isppol)=buffer(jdim+1:jdim+dimlmn(iatom,irhoij)) 3040 jdim=jdim+dimlmn(iatom,irhoij) 3041 end do 3042 end do 3043 end do 3044 3045 LIBPAW_DEALLOCATE(buffer) 3046 LIBPAW_DEALLOCATE(dimlmn) 3047 3048 end subroutine pawrhoij_mpisum_unpacked_2D
m_pawrhoij/pawrhoij_nullify [ Functions ]
[ Top ] [ m_pawrhoij ] [ Functions ]
NAME
pawrhoij_nullify
FUNCTION
Nullify (initialize to null) a pawrhoij datastructure
SIDE EFFECTS
pawrhoij(:)<type(pawrhoij_type)>= rhoij datastructure
SOURCE
448 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 nrhoij=size(pawrhoij) 465 466 if (nrhoij>0) then 467 do irhoij=1,nrhoij 468 pawrhoij(irhoij)%cplex_rhoij=1 469 pawrhoij(irhoij)%qphase=1 470 pawrhoij(irhoij)%nrhoijsel=0 471 pawrhoij(irhoij)%ngrhoij=0 472 pawrhoij(irhoij)%lmnmix_sz=0 473 pawrhoij(irhoij)%use_rhoij_=0 474 pawrhoij(irhoij)%use_rhoijp=0 475 pawrhoij(irhoij)%use_rhoijres=0 476 end do 477 end if 478 479 end subroutine pawrhoij_nullify
m_pawrhoij/pawrhoij_redistribute [ Functions ]
[ Top ] [ m_pawrhoij ] [ Functions ]
NAME
pawrhoij_redistribute
FUNCTION
Redistribute an array of pawrhoij datastructures Input pawrhoij is given on a MPI communicator Output pawrhoij is redistributed on another MPI communicator
INPUTS
mpi_comm_in= input MPI (atom) communicator
mpi_comm_out= output MPI (atom) communicator
mpi_atmtab_in= --optional-- indexes of the input pawrhoij treated by current proc
if not present, will be calculated in the present routine
mpi_atmtab_out= --optional-- indexes of the output pawrhoij treated by current proc
if not present, will be calculated in the present routine
natom= --optional-- total number of atoms
----- Optional arguments used only for asynchronous communications -----
RecvAtomProc(:)= rank of processor from which I expect atom (in mpi_comm_in)
RecvAtomList(:)= indexes of atoms to be received by me
RecvAtomList(irecv) are the atoms I expect from RecvAtomProc(irecv)
SendAtomProc(:)= ranks of process destination of atom (in mpi_comm_in)
SendAtomList(:)= indexes of atoms to be sent by me
SendAtomList(isend) are the atoms sent to SendAtomProc(isend)
OUTPUT
[pawrhoij_out(:)]<type(pawrhoij_type)>= --optional-- if present, the redistributed datastructure does not replace the input one but is delivered in pawrhoij_out if not present, input and output datastructure are the same.
SIDE EFFECTS
pawrhoij(:)<type(pawrhoij_type)>= input (and eventually output) pawrhoij datastructures
SOURCE
1925 subroutine pawrhoij_redistribute(pawrhoij,mpi_comm_in,mpi_comm_out,& 1926 & natom,mpi_atmtab_in,mpi_atmtab_out,pawrhoij_out,& 1927 & 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 integer,pointer :: my_atmtab_in(:),my_atmtab_out(:) 1949 integer,allocatable :: atmtab_send(:),atm_indx_in(:),atm_indx_out(:),buf_int1(:),From(:),request(:) 1950 integer,allocatable,target :: buf_int(:) 1951 integer,pointer:: buf_ints(:) 1952 logical,allocatable :: msg_pick(:) 1953 real(dp),allocatable :: buf_dp1(:) 1954 real(dp),allocatable,target :: buf_dp(:) 1955 real(dp),pointer :: buf_dps(:) 1956 type(coeffi1_type),target,allocatable :: tab_buf_int(:),tab_buf_atom(:) 1957 type(coeff1_type),target,allocatable :: tab_buf_dp(:) 1958 type(pawrhoij_type),pointer :: pawrhoij_out1(:) 1959 type(pawrhoij_type),allocatable :: pawrhoij_all(:) 1960 1961 ! ************************************************************************* 1962 1963 !@pawrhoij_type 1964 1965 in_place=(.not.present(pawrhoij_out)) 1966 my_natom_in=size(pawrhoij) 1967 1968 !If not "in_place", destroy the output datastructure 1969 if (.not.in_place) then 1970 if (associated(pawrhoij_out)) then 1971 call pawrhoij_free(pawrhoij_out) 1972 LIBPAW_DATATYPE_DEALLOCATE(pawrhoij_out) 1973 end if 1974 end if 1975 1976 !Special sequential case 1977 if (mpi_comm_in==xmpi_comm_self.and.mpi_comm_out==xmpi_comm_self) then 1978 if ((.not.in_place).and.(my_natom_in>0)) then 1979 LIBPAW_DATATYPE_ALLOCATE(pawrhoij_out,(my_natom_in)) 1980 call pawrhoij_nullify(pawrhoij_out) 1981 call pawrhoij_copy(pawrhoij,pawrhoij_out,& 1982 & keep_cplex=.false.,keep_qphase=.false.,keep_itypat=.false.,keep_nspden=.false.) 1983 end if 1984 return 1985 end if 1986 1987 !Get total natom 1988 if (present(natom)) then 1989 natom_tot=natom 1990 else 1991 natom_tot=my_natom_in 1992 call xmpi_sum(natom_tot,mpi_comm_in,ierr) 1993 end if 1994 1995 !Select input distribution 1996 if (present(mpi_atmtab_in)) then 1997 my_atmtab_in => mpi_atmtab_in 1998 my_atmtab_in_allocated=.false. 1999 else 2000 call get_my_atmtab(mpi_comm_in,my_atmtab_in,my_atmtab_in_allocated,& 2001 & paral_atom,natom_tot,my_natom_in) 2002 end if 2003 2004 !Select output distribution 2005 if (present(mpi_atmtab_out)) then 2006 my_natom_out=size(mpi_atmtab_out) 2007 my_atmtab_out => mpi_atmtab_out 2008 my_atmtab_out_allocated=.false. 2009 else 2010 call get_my_atmtab(mpi_comm_out,my_atmtab_out,my_atmtab_out_allocated,& 2011 & paral_atom,natom_tot) 2012 end if 2013 2014 !Select algo according to optional input arguments 2015 algo_option=1 2016 if (present(SendAtomProc).and.present(SendAtomList).and.& 2017 & 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 LIBPAW_DATATYPE_ALLOCATE(pawrhoij_all,(natom_tot)) 2025 call pawrhoij_nullify(pawrhoij_all) 2026 call pawrhoij_copy(pawrhoij,pawrhoij_all,comm_atom=mpi_comm_in,mpi_atmtab=my_atmtab_in,& 2027 & keep_cplex=.false.,keep_qphase=.false.,keep_itypat=.false.,keep_nspden=.false.) 2028 if (in_place) then 2029 call pawrhoij_free(pawrhoij) 2030 LIBPAW_DATATYPE_DEALLOCATE(pawrhoij) 2031 LIBPAW_DATATYPE_ALLOCATE(pawrhoij,(my_natom_out)) 2032 call pawrhoij_nullify(pawrhoij) 2033 call pawrhoij_copy(pawrhoij_all,pawrhoij,comm_atom=mpi_comm_out,mpi_atmtab=my_atmtab_out,& 2034 & keep_cplex=.false.,keep_qphase=.false.,keep_itypat=.false.,keep_nspden=.false.) 2035 else 2036 LIBPAW_DATATYPE_ALLOCATE(pawrhoij_out,(my_natom_out)) 2037 call pawrhoij_nullify(pawrhoij_out) 2038 call pawrhoij_copy(pawrhoij_all,pawrhoij_out,comm_atom=mpi_comm_out,mpi_atmtab=my_atmtab_out,& 2039 & keep_cplex=.false.,keep_qphase=.false.,keep_itypat=.false.,keep_nspden=.false.) 2040 end if 2041 call pawrhoij_free(pawrhoij_all) 2042 LIBPAW_DATATYPE_DEALLOCATE(pawrhoij_all) 2043 2044 2045 !Asynchronous algorithm (asynchronous communications) 2046 !--------------------------------------------------------- 2047 else if (algo_option==2) then 2048 2049 nbsend=size(SendAtomProc) ; nbrecv=size(RecvAtomProc) 2050 2051 if (in_place) then 2052 if ( my_natom_out > 0 ) then 2053 LIBPAW_DATATYPE_ALLOCATE(pawrhoij_out1,(my_natom_out)) 2054 call pawrhoij_nullify(pawrhoij_out1) 2055 else 2056 LIBPAW_DATATYPE_ALLOCATE(pawrhoij_out1,(0)) 2057 end if 2058 else 2059 LIBPAW_DATATYPE_ALLOCATE(pawrhoij_out,(my_natom_out)) 2060 call pawrhoij_nullify(pawrhoij_out) 2061 pawrhoij_out1=>pawrhoij_out 2062 end if 2063 2064 nproc_in=xmpi_comm_size(mpi_comm_in) 2065 nproc_out=xmpi_comm_size(mpi_comm_out) 2066 if (nproc_in<=nproc_out) mpi_comm_exch=mpi_comm_out 2067 if (nproc_in>nproc_out) mpi_comm_exch=mpi_comm_in 2068 me_exch=xmpi_comm_rank(mpi_comm_exch) 2069 2070 ! Dimension put to the maximum to send 2071 LIBPAW_ALLOCATE(atmtab_send,(nbsend)) 2072 LIBPAW_ALLOCATE(atm_indx_in,(natom_tot)) 2073 atm_indx_in=-1 2074 do iatom=1,my_natom_in 2075 atm_indx_in(my_atmtab_in(iatom))=iatom 2076 end do 2077 LIBPAW_ALLOCATE(atm_indx_out,(natom_tot)) 2078 atm_indx_out=-1 2079 do iatom=1,my_natom_out 2080 atm_indx_out(my_atmtab_out(iatom))=iatom 2081 end do 2082 2083 LIBPAW_DATATYPE_ALLOCATE(tab_buf_int,(nbsend)) 2084 LIBPAW_DATATYPE_ALLOCATE(tab_buf_dp,(nbsend)) 2085 LIBPAW_DATATYPE_ALLOCATE(tab_buf_atom,(nbsend)) 2086 LIBPAW_ALLOCATE(request,(3*nbsend)) 2087 2088 ! A send buffer in an asynchrone communication couldn't be deallocate before it has been receive 2089 nbsent=0 ; ireq=0 ; iisend=0 ; nbsendreq=0 ; nb_msg=0 2090 do iisend=1,nbsend 2091 iproc_rcv=SendAtomProc(iisend) 2092 next=-1 2093 if (iisend < nbsend) next=SendAtomProc(iisend+1) 2094 if (iproc_rcv /= me_exch) then 2095 nbsent=nbsent+1 2096 atmtab_send(nbsent)=SendAtomList(iisend) ! we groups the atoms sends to the same process 2097 if (iproc_rcv /= next) then 2098 if (nbsent > 0) then 2099 ! Check if message has been yet prepared 2100 message_yet_prepared=.false. 2101 do imsg=1,nb_msg 2102 if (size(tab_buf_atom(imsg)%value) /= nbsent) then 2103 cycle 2104 else 2105 do imsg1=1,nbsent 2106 if (tab_buf_atom(imsg)%value(imsg1)/=atmtab_send(imsg1)) exit 2107 message_yet_prepared=.true. 2108 imsg_current=imsg 2109 end do 2110 end if 2111 end do 2112 ! Create the message 2113 if (.not.message_yet_prepared) then 2114 nb_msg=nb_msg+1 2115 call pawrhoij_isendreceive_fillbuffer( & 2116 & pawrhoij,atmtab_send,atm_indx_in,nbsent,buf_int,nb_int,buf_dp,nb_dp) 2117 LIBPAW_ALLOCATE(tab_buf_int(nb_msg)%value,(nb_int)) 2118 LIBPAW_ALLOCATE(tab_buf_dp(nb_msg)%value,(nb_dp)) 2119 tab_buf_int(nb_msg)%value(1:nb_int)=buf_int(1:nb_int) 2120 tab_buf_dp(nb_msg)%value(1:nb_dp)=buf_dp(1:nb_dp) 2121 LIBPAW_DEALLOCATE(buf_int) 2122 LIBPAW_DEALLOCATE(buf_dp) 2123 LIBPAW_ALLOCATE(tab_buf_atom(nb_msg)%value, (nbsent)) 2124 tab_buf_atom(nb_msg)%value(1:nbsent)=atmtab_send(1:nbsent) 2125 imsg_current=nb_msg 2126 end if 2127 ! Communicate 2128 buf_size(1)=size(tab_buf_int(imsg_current)%value) 2129 buf_size(2)=size(tab_buf_dp(imsg_current)%value) 2130 buf_size(3)=nbsent 2131 buf_ints=>tab_buf_int(imsg_current)%value 2132 buf_dps=>tab_buf_dp(imsg_current)%value 2133 my_tag=100 2134 ireq=ireq+1 2135 call xmpi_isend(buf_size,iproc_rcv,my_tag,mpi_comm_exch,request(ireq),ierr) 2136 my_tag=101 2137 ireq=ireq+1 2138 call xmpi_isend(buf_ints,iproc_rcv,my_tag,mpi_comm_exch,request(ireq),ierr) 2139 my_tag=102 2140 ireq=ireq+1 2141 call xmpi_isend(buf_dps,iproc_rcv,my_tag,mpi_comm_exch,request(ireq),ierr) 2142 nbsendreq=ireq 2143 nbsent=0 2144 end if 2145 end if 2146 else ! Just a renumbering, not a sending 2147 iat_in=atm_indx_in(SendAtomList(iisend)) 2148 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 & keep_cplex=.false.,keep_qphase=.false.,keep_itypat=.false.,keep_nspden=.false.) 2151 nbsent=0 2152 end if 2153 end do 2154 2155 LIBPAW_ALLOCATE(From,(nbrecv)) 2156 From(:)=-1 ; nbrecvmsg=0 2157 do iircv=1,nbrecv 2158 iproc_send=RecvAtomProc(iircv) !receive from (RcvAtomProc is sorted by growing process) 2159 next=-1 2160 if (iircv < nbrecv) next=RecvAtomProc(iircv+1) 2161 if (iproc_send /= me_exch .and. iproc_send/=next) then 2162 nbrecvmsg=nbrecvmsg+1 2163 From(nbrecvmsg)=iproc_send 2164 end if 2165 end do 2166 2167 LIBPAW_ALLOCATE(msg_pick,(nbrecvmsg)) 2168 msg_pick=.false. 2169 nbmsg_incoming=nbrecvmsg 2170 do while (nbmsg_incoming > 0) 2171 do i1=1,nbrecvmsg 2172 if (.not.msg_pick(i1)) then 2173 iproc_send=From(i1) 2174 flag=.false. 2175 my_tag=100 2176 call xmpi_iprobe(iproc_send,my_tag,mpi_comm_exch,flag,ierr) 2177 if (flag) then 2178 msg_pick(i1)=.true. 2179 call xmpi_irecv(buf_size,iproc_send,my_tag,mpi_comm_exch,request1(1),ierr) 2180 call xmpi_wait(request1(1),ierr) 2181 nb_int=buf_size(1) 2182 nb_dp=buf_size(2) 2183 npawrhoij_sent=buf_size(3) 2184 LIBPAW_ALLOCATE(buf_int1,(nb_int)) 2185 LIBPAW_ALLOCATE(buf_dp1,(nb_dp)) 2186 my_tag=101 2187 call xmpi_irecv(buf_int1,iproc_send,my_tag,mpi_comm_exch,request1(2),ierr) 2188 my_tag=102 2189 call xmpi_irecv(buf_dp1,iproc_send,my_tag,mpi_comm_exch,request1(3),ierr) 2190 call xmpi_waitall(request1(2:3),ierr) 2191 call pawrhoij_isendreceive_getbuffer(pawrhoij_out1,npawrhoij_sent,atm_indx_out,buf_int1,buf_dp1) 2192 nbmsg_incoming=nbmsg_incoming-1 2193 LIBPAW_DEALLOCATE(buf_int1) 2194 LIBPAW_DEALLOCATE(buf_dp1) 2195 end if 2196 end if 2197 end do 2198 end do 2199 LIBPAW_DEALLOCATE(msg_pick) 2200 2201 if (in_place) then 2202 call pawrhoij_free(pawrhoij) 2203 LIBPAW_DATATYPE_DEALLOCATE(pawrhoij) 2204 LIBPAW_DATATYPE_ALLOCATE(pawrhoij,(my_natom_out)) 2205 call pawrhoij_nullify(pawrhoij) 2206 call pawrhoij_copy(pawrhoij_out1,pawrhoij, & 2207 & keep_cplex=.false.,keep_qphase=.false.,keep_itypat=.false.,keep_nspden=.false.) 2208 call pawrhoij_free(pawrhoij_out1) 2209 LIBPAW_DATATYPE_DEALLOCATE(pawrhoij_out1) 2210 end if 2211 2212 ! Wait for deallocating arrays that all sending operations has been realized 2213 if (nbsendreq > 0) then 2214 call xmpi_waitall(request(1:nbsendreq),ierr) 2215 end if 2216 2217 ! Deallocate buffers 2218 do i1=1,nb_msg 2219 LIBPAW_DEALLOCATE(tab_buf_int(i1)%value) 2220 LIBPAW_DEALLOCATE(tab_buf_dp(i1)%value) 2221 LIBPAW_DEALLOCATE(tab_buf_atom(i1)%value) 2222 end do 2223 LIBPAW_DATATYPE_DEALLOCATE(tab_buf_int) 2224 LIBPAW_DATATYPE_DEALLOCATE(tab_buf_dp) 2225 LIBPAW_DATATYPE_DEALLOCATE(tab_buf_atom) 2226 LIBPAW_DEALLOCATE(From) 2227 LIBPAW_DEALLOCATE(request) 2228 LIBPAW_DEALLOCATE(atmtab_send) 2229 LIBPAW_DEALLOCATE(atm_indx_in) 2230 LIBPAW_DEALLOCATE(atm_indx_out) 2231 2232 end if !algo_option 2233 2234 !Eventually release temporary pointers 2235 call free_my_atmtab(my_atmtab_in,my_atmtab_in_allocated) 2236 call free_my_atmtab(my_atmtab_out,my_atmtab_out_allocated) 2237 2238 end subroutine pawrhoij_redistribute
m_pawrhoij/pawrhoij_symrhoij [ Functions ]
[ Top ] [ m_pawrhoij ] [ Functions ]
NAME
pawrhoij_symrhoij
FUNCTION
Symmetrize rhoij quantities (augmentation occupancies) and/or gradients Compute also rhoij residuals (new-old values of rhoij and gradients)
INPUTS
choice=select then type of rhoij gradients to symmetrize.
choice=1 => no gradient
choice=2 => gradient with respect to atomic position(s)
=3 => a gradient with respect to strain(s)
=4 => 2nd gradient with respect to atomic position(s)
=23=> a gradient with respect to atm. pos. and strain(s)
=24=> 1st and 2nd gradient with respect to atomic position(s)
gprimd(3,3)=dimensional primitive translations for reciprocal space(bohr^-1).
indsym(4,nsym,natom)=indirect indexing array for atom labels
ipert=index of perturbation if pawrhoij is a pertubed rhoij
no meaning for ground-state calculations (should be 0)
[mpi_atmtab(:)]=--optional-- indexes of the atoms treated by current proc
[comm_atom]=--optional-- MPI communicator over atoms
natom=number of atoms in cell
nsym=number of symmetry elements in space group
ntypat=number of types of atoms in unit cell.
optrhoij= 1 if rhoij quantities have to be symmetrized
pawrhoij_unsym(:) <type(pawrhoij_type)>= datastructure containing PAW rhoij occupancies
(and related data) non symmetrized in an unpacked storage (pawrhoij_unsym%rhoij_)
Contains eventually unsymmetrized rhoij gradients (grhoij)
pawang <type(pawang_type)>=angular mesh discretization and related data
pawprtvol=control print volume and debugging output for PAW
Note: if pawprtvol=-10001, nothing is printed out
pawtab(ntypat) <type(pawtab_type)>=paw tabulated starting data
[qphon(3)]=--optional-- (RF calculations only) - wavevector of the phonon
rprimd(3,3)=real space primitive translations.
symafm(nsym)=(anti)ferromagnetic part of symmetry operations
symrec(3,3,nsym)=symmetries of group in terms of operations on
reciprocal space primitive translations
typat(natom)=type for each atom
[use_zeromag]=--optional-- .TRUE. if rhoij "magnetization" is enforced to be zero
Applies only when nspden_rhoij=4 (note: only the real part is set to zero)
OUTPUT
SIDE EFFECTS
pawrhoij(natom) <type(pawrhoij_type)>= datastructure containing PAW rhoij occupancies (and related data) SYMMETRIZED in a packed storage (pawrhoij%rhoijp). if (optrhoij==1) pawrhoij(natom)%nrhoijsel=number of non-zero values of rhoij pawrhoij(natom)%rhoijp(:,:)=symetrized paw rhoij quantities in PACKED STORAGE (only non-zero values) pawrhoij(natom)%rhoijres(:,:)=paw rhoij quantities residuals (new values - old values) pawrhoij(natom)%rhoijselect(:)=select the non-zero values of rhoij!! if (pawrhoij(:)%ngrhoij>0) (equivalent to choice>1) pawrhoij(natom)%grhoij(:,:)=symetrized gradients of rhoij
NOTES
pawrhoij and pawrhoij_unsym can be identical (refer to the same pawrhoij datastructure). They should be different only if pawrhoij is distributed over atomic sites (in that case pawrhoij_unsym should not be distributed over atomic sites).
SOURCE
3491 subroutine pawrhoij_symrhoij(pawrhoij,pawrhoij_unsym,choice,gprimd,indsym,ipert,natom,nsym,& 3492 & ntypat,optrhoij,pawang,pawprtvol,pawtab,rprimd,symafm,symrec,typat, & 3493 & 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 logical,parameter :: afm_noncoll=.true. ! TRUE if antiferro symmetries are used with non-collinear magnetism 3521 logical :: use_zeromag_ 3522 real(dp) :: arg,factafm,ro,syma,zarot2 3523 logical :: antiferro,has_qphase,my_atmtab_allocated,noncoll 3524 logical :: paral_atom,paral_atom_unsym,use_afm,use_res 3525 character(len=8) :: pertstrg,wrt_mode 3526 character(len=500) :: msg 3527 !arrays 3528 integer,parameter :: alpha(6)=(/1,2,3,3,3,2/),beta(6)=(/1,2,3,2,1,1/) 3529 integer :: nsym_used(2) 3530 integer, pointer :: indlmn(:,:) 3531 integer,pointer :: my_atmtab(:) 3532 real(dp) :: fact(2),factsym(2),phase(2),rhoijc(2),rotmag(2,3,2),rotrho(2,2,2) 3533 real(dp) :: summag(2,3,2),sumrho(2,2,2),sum1(2),work1(2,3,3),xsym(3) 3534 real(dp),allocatable :: rotgr(:,:,:,:),rotmaggr(:,:,:,:),sumgr(:,:,:),summaggr(:,:,:,:) 3535 real(dp),allocatable :: symrec_cart(:,:,:) 3536 real(dp),pointer :: grad(:,:,:) 3537 type(coeff3_type),target,allocatable :: tmp_grhoij(:) 3538 type(pawrhoij_type),pointer :: pawrhoij_unsym_all(:) 3539 3540 ! ********************************************************************* 3541 3542 !Sizes of pawrhoij datastructures 3543 nrhoij=size(pawrhoij) 3544 nrhoij_unsym=size(pawrhoij_unsym) 3545 3546 !Set up parallelism over atoms 3547 paral_atom=(present(comm_atom).and.(nrhoij/=natom)) 3548 paral_atom_unsym=(present(comm_atom).and.(nrhoij_unsym/=natom)) 3549 nullify(my_atmtab);if (present(mpi_atmtab)) my_atmtab => mpi_atmtab 3550 my_comm_atom=xmpi_comm_self;if (present(comm_atom)) my_comm_atom=comm_atom 3551 call get_my_atmtab(my_comm_atom,my_atmtab,my_atmtab_allocated,paral_atom,natom) 3552 3553 !Test: consistency between choice and ngrhoij 3554 ngrhoij=0 3555 if (nrhoij>0) then 3556 ngrhoij=pawrhoij(1)%ngrhoij 3557 if(choice==2) ngrhoij=min(3,pawrhoij(1)%ngrhoij) 3558 if(choice==3.or.choice==4) ngrhoij=min(6,pawrhoij(1)%ngrhoij) 3559 if(choice==23.or.choice==24) ngrhoij=min(9,pawrhoij(1)%ngrhoij) 3560 if ((choice==1.and.ngrhoij/=0).or.(choice==2.and.ngrhoij/=3).or. & 3561 & (choice==3.and.ngrhoij/=6).or.(choice==23.and.ngrhoij/=9).or. & 3562 & (choice==4.and.ngrhoij/=6).or.(choice==24.and.ngrhoij/=9) ) then 3563 msg='Inconsistency between variables choice and ngrhoij !' 3564 LIBPAW_BUG(msg) 3565 end if 3566 end if 3567 3568 !Antiferro case ? 3569 antiferro=.false.;if (nrhoij>0) antiferro=(pawrhoij(1)%nspden==2.and.pawrhoij(1)%nsppol==1) 3570 !Non-collinear case 3571 noncoll=.false.;if (nrhoij>0) noncoll=(pawrhoij(1)%nspden==4) 3572 !Do we use antiferro symmetries ? 3573 use_afm=((antiferro).or.(noncoll.and.afm_noncoll)) 3574 !Do we impose zero magnetization? 3575 use_zeromag_=.false. ; if (present(use_zeromag)) use_zeromag_=use_zeromag 3576 3577 ! Does not symmetrize imaginary part for GS calculations 3578 cplex_eff=1 3579 if (nrhoij>0.and.(ipert>0.or.antiferro.or.noncoll)) cplex_eff=pawrhoij(1)%cplex_rhoij 3580 3581 !Do we have a phase due to q-vector? 3582 has_qphase=.false. 3583 if (nrhoij>0) then 3584 has_qphase=(pawrhoij(1)%qphase==2) 3585 if (present(qphon)) then 3586 if (any(abs(qphon(1:3))>tol8).and.(.not.has_qphase)) then 3587 msg='Should have qphase=2 for a non-zero q!' 3588 LIBPAW_BUG(msg) 3589 end if 3590 end if 3591 end if 3592 3593 !Printing of unsymetrized Rhoij 3594 if (nrhoij>0.and.optrhoij==1.and.pawprtvol/=-10001) then 3595 wrt_mode='COLL';if (paral_atom) wrt_mode='PERS' 3596 pertstrg="RHOIJ";if (ipert>0) pertstrg="RHOIJ(1)" 3597 natinc=1;if(nrhoij>1.and.pawprtvol>=0) natinc=nrhoij-1 3598 write(msg, '(7a)') ch10," PAW TEST:",ch10,& 3599 & ' ========= Values of ',trim(pertstrg),' before symetrization =========',ch10 3600 call wrtout(std_out,msg,wrt_mode) 3601 do iatm=1,nrhoij,natinc 3602 iatom=iatm; if (paral_atom) iatom=my_atmtab(iatm) 3603 if (nrhoij==1.and.ipert>0.and.ipert<=natom) iatom=ipert 3604 call pawrhoij_print_rhoij(pawrhoij_unsym(iatm)%rhoij_,pawrhoij_unsym(iatm)%cplex_rhoij,& 3605 & pawrhoij_unsym(iatm)%qphase,iatom,natom,& 3606 & unit=std_out,opt_prtvol=pawprtvol,mode_paral=wrt_mode) 3607 end do 3608 call wrtout(std_out,"",wrt_mode) 3609 end if 3610 3611 !Symetrization occurs only when nsym>1 3612 if (nsym>1) then 3613 3614 ! Symetrization of gradients not compatible with nspden=4 3615 if (nrhoij>0) then 3616 if (choice>2.and.pawrhoij(1)%nspden==4) then 3617 msg='For the time being, choice>2 is not compatible with nspden=4 !' 3618 LIBPAW_BUG(msg) 3619 end if 3620 end if 3621 3622 ! Symetry matrixes must be in memory 3623 if (pawang%nsym==0) then 3624 msg='pawang%zarot must be allocated !' 3625 LIBPAW_BUG(msg) 3626 end if 3627 3628 if (has_qphase.and.choice>1) then 3629 msg='choice>1 not compatible with q-phase !' 3630 LIBPAW_BUG(msg) 3631 end if 3632 3633 ! Several inits/allocations 3634 if (noncoll) then 3635 LIBPAW_ALLOCATE(symrec_cart,(3,3,nsym)) 3636 do irot=1,nsym 3637 symrec_cart(:,:,irot)=symrhoij_symcart(gprimd,rprimd,symrec(:,:,irot)) 3638 end do 3639 end if 3640 ishift2=0;ishift3=0;ishift4=0 3641 if (choice>1) then 3642 iafm=merge(2,1,antiferro) 3643 qphase=merge(2,1,has_qphase) 3644 LIBPAW_ALLOCATE(sumgr,(cplex_eff,ngrhoij,qphase)) 3645 LIBPAW_ALLOCATE(rotgr,(cplex_eff,ngrhoij,iafm,qphase)) 3646 if (noncoll) then 3647 LIBPAW_ALLOCATE(summaggr,(cplex_eff,ngrhoij,3,qphase)) 3648 LIBPAW_ALLOCATE(rotmaggr,(cplex_eff,ngrhoij,3,qphase)) 3649 end if 3650 if (choice==23) ishift2=6 3651 if (choice==24) ishift4=3 3652 if (.not.paral_atom_unsym) then 3653 ! Have to make a temporary copy of grhoij 3654 LIBPAW_DATATYPE_ALLOCATE(tmp_grhoij,(nrhoij)) 3655 do iatm=1,nrhoij 3656 sz1=pawrhoij_unsym(iatm)%cplex_rhoij*pawrhoij_unsym(iatm)%qphase & 3657 & *pawrhoij_unsym(iatm)%lmn2_size 3658 sz2=pawrhoij_unsym(iatm)%nspden 3659 LIBPAW_ALLOCATE(tmp_grhoij(iatm)%value,(ngrhoij,sz1,sz2)) 3660 tmp_grhoij(iatm)%value(1:ngrhoij,1:sz1,1:sz2)= & 3661 & pawrhoij_unsym(iatm)%grhoij(1:ngrhoij,1:sz1,1:sz2) 3662 end do 3663 end if 3664 end if 3665 3666 ! In case of paw_rhoij_unsym distributed over atomic sites, gather it 3667 if (paral_atom_unsym) then 3668 LIBPAW_DATATYPE_ALLOCATE(pawrhoij_unsym_all,(natom)) 3669 call pawrhoij_nullify(pawrhoij_unsym_all) 3670 call pawrhoij_gather(pawrhoij_unsym,pawrhoij_unsym_all,-1,my_comm_atom,& 3671 & with_lmnmix=.false.,with_rhoijp=.false.,& 3672 & with_rhoijres=.false.,with_grhoij=(choice>1)) 3673 nrhoij1=natom 3674 else 3675 pawrhoij_unsym_all=>pawrhoij_unsym 3676 nrhoij1=nrhoij_unsym 3677 end if 3678 3679 3680 ! Loops over atoms and spin components 3681 ! ------------------------------------ 3682 do iatm=1,nrhoij 3683 iatom=iatm;if (paral_atom) iatom=my_atmtab(iatm) 3684 if (nrhoij==1.and.ipert>0.and.ipert<=natom.and.(.not.paral_atom)) iatom=ipert 3685 itypat=typat(iatom) 3686 lmn_size=pawrhoij(iatm)%lmn_size 3687 lmn2_size=pawrhoij(iatm)%lmn2_size 3688 qphase=pawrhoij(iatm)%qphase 3689 cplex_rhoij=pawrhoij(iatm)%cplex_rhoij 3690 cplex_eff=1;if (ipert>0.or.antiferro.or.noncoll) cplex_eff=cplex_rhoij 3691 use_res=(pawrhoij(iatm)%use_rhoijres>0) 3692 indlmn => pawtab(itypat)%indlmn 3693 3694 nselect=0 3695 do ispden=1,pawrhoij(iatm)%nsppol 3696 jspden=min(3-ispden,pawrhoij(iatm)%nsppol) 3697 3698 ! Store old -rhoij in residual 3699 if (optrhoij==1.and.use_res) then 3700 pawrhoij(iatm)%rhoijres(:,ispden)=zero 3701 if (noncoll) pawrhoij(iatm)%rhoijres(:,2:4)=zero 3702 if (antiferro) pawrhoij(iatm)%rhoijres(:,2)=zero 3703 do iq=1,qphase 3704 iq0=(iq-1)*cplex_rhoij*lmn2_size 3705 if (cplex_rhoij==1) then 3706 do irhoij=1,pawrhoij(iatm)%nrhoijsel 3707 klmn1=iq0+pawrhoij(iatm)%rhoijselect(irhoij);jrhoij=iq0+irhoij 3708 pawrhoij(iatm)%rhoijres(klmn1,ispden)=-pawrhoij(iatm)%rhoijp(jrhoij,ispden) 3709 end do 3710 else 3711 do irhoij=1,pawrhoij(iatm)%nrhoijsel 3712 klmn1=iq0+2*pawrhoij(iatm)%rhoijselect(irhoij);jrhoij=iq0+2*irhoij 3713 pawrhoij(iatm)%rhoijres(klmn1-1:klmn1,ispden)=-pawrhoij(iatm)%rhoijp(jrhoij-1:jrhoij,ispden) 3714 end do 3715 end if 3716 if (noncoll) then 3717 if (cplex_rhoij==1) then 3718 do mu=2,4 3719 do irhoij=1,pawrhoij(iatm)%nrhoijsel 3720 klmn1=iq0+pawrhoij(iatm)%rhoijselect(irhoij);jrhoij=iq0+irhoij 3721 pawrhoij(iatm)%rhoijres(klmn1,mu)=-pawrhoij(iatm)%rhoijp(jrhoij,mu) 3722 end do 3723 end do 3724 else 3725 do mu=2,4 3726 do irhoij=1,pawrhoij(iatm)%nrhoijsel 3727 klmn1=iq0+2*pawrhoij(iatm)%rhoijselect(irhoij);jrhoij=iq0+2*irhoij 3728 pawrhoij(iatm)%rhoijres(klmn1-1:klmn1,mu)=-pawrhoij(iatm)%rhoijp(jrhoij-1:jrhoij,mu) 3729 end do 3730 end do 3731 end if 3732 end if 3733 if (antiferro) then 3734 if (cplex_rhoij==1) then 3735 do irhoij=1,pawrhoij(iatm)%nrhoijsel 3736 klmn1=iq0+pawrhoij(iatm)%rhoijselect(irhoij);jrhoij=iq0+irhoij 3737 pawrhoij(iatm)%rhoijres(klmn1,2)=-pawrhoij(iatm)%rhoijp(jrhoij,2) 3738 end do 3739 else 3740 do irhoij=1,pawrhoij(iatm)%nrhoijsel 3741 klmn1=iq0+2*pawrhoij(iatm)%rhoijselect(irhoij);jrhoij=iq0+2*irhoij 3742 pawrhoij(iatm)%rhoijres(klmn1-1:klmn1,2)=-pawrhoij(iatm)%rhoijp(jrhoij-1:jrhoij,2) 3743 end do 3744 end if 3745 end if 3746 end do 3747 end if 3748 3749 3750 ! Loops over (il,im) and (jl,jm) 3751 ! ------------------------------ 3752 jl0=-1;jln0=-1;indexj=1 3753 do jlmn=1,lmn_size 3754 jl=indlmn(1,jlmn) 3755 jlpm=1+jl+indlmn(2,jlmn) 3756 jln=indlmn(5,jlmn) 3757 if (jln/=jln0) indexj=indexj+2*jl0+1 3758 j0lmn=jlmn*(jlmn-1)/2 3759 il0=-1;iln0=-1;indexi=1 3760 do ilmn=1,jlmn 3761 il=indlmn(1,ilmn) 3762 ilpm=1+il+indlmn(2,ilmn) 3763 iln=indlmn(5,ilmn) 3764 if (iln/=iln0) indexi=indexi+2*il0+1 3765 klmn=j0lmn+ilmn 3766 klmn1=merge(klmn,2*klmn-1,cplex_rhoij==1) 3767 3768 nsym_used(:)=0 3769 if (optrhoij==1) rotrho=zero 3770 if (optrhoij==1.and.noncoll) rotmag=zero 3771 if (choice>1) rotgr=zero 3772 if (choice>1.and.noncoll) rotmaggr=zero 3773 3774 3775 ! Loop over symmetries 3776 ! -------------------- 3777 do irot=1,nsym 3778 3779 if ((symafm(irot)/=1).and.(.not.use_afm)) cycle 3780 kspden=ispden;if (symafm(irot)==-1) kspden=jspden 3781 iafm=1;if ((antiferro).and.(symafm(irot)==-1)) iafm=2 3782 factafm=dble(symafm(irot)) 3783 3784 nsym_used(iafm)=nsym_used(iafm)+1 3785 at_indx=min(indsym(4,irot,iatom),nrhoij1) 3786 3787 if (has_qphase) then 3788 arg=two_pi*(qphon(1)*indsym(1,irot,iatom)+qphon(2)*indsym(2,irot,iatom) & 3789 & +qphon(3)*indsym(3,irot,iatom)) 3790 phase(1)=cos(arg);phase(2)=sin(arg) 3791 end if 3792 3793 if (optrhoij==1) sumrho=zero 3794 if (optrhoij==1.and.noncoll) summag=zero 3795 if (choice>1) sumgr=zero 3796 if (choice>1.and.noncoll) summaggr=zero 3797 3798 if (choice>1) then 3799 if (paral_atom_unsym) then 3800 grad => pawrhoij_unsym_all(at_indx)%grhoij 3801 else 3802 grad => tmp_grhoij(at_indx)%value 3803 end if 3804 end if 3805 3806 3807 ! Accumulate values over (mi,mj) 3808 ! ------------------------------ 3809 do mj=1,2*jl+1 3810 indexjj=indexj+mj;indexjj0=indexjj*(indexjj-1)/2 3811 do mi=1,2*il+1 3812 factsym(:)=one 3813 indexii=indexi+mi 3814 if (indexii<=indexjj) then 3815 indexk=indexjj0+indexii 3816 factsym(2)=one 3817 else 3818 indexk=indexii*(indexii-1)/2+indexjj 3819 factsym(2)=-one 3820 end if 3821 indexkc=cplex_rhoij*(indexk-1) 3822 indexkc_q=indexkc+cplex_rhoij*lmn2_size 3823 3824 ! Be careful: use here R_rel^-1 in term of spherical harmonics 3825 ! which is tR_rec in term of spherical harmonics 3826 ! so, use transpose[zarot] 3827 zarot2=pawang%zarot(mi,ilpm,il+1,irot)*pawang%zarot(mj,jlpm,jl+1,irot) 3828 ! zarot2=pawang%zarot(ilpm,mi,il+1,irot)*pawang%zarot(jlpm,mj,jl+1,irot) 3829 3830 ! Rotate rhoij 3831 if (optrhoij==1) then 3832 fact(1)=factsym(1);fact(2)=factsym(2)*factafm !????? What? MT 3833 sumrho(1:cplex_eff,iafm,1)=sumrho(1:cplex_eff,iafm,1) & 3834 & +fact(1:cplex_eff)*zarot2 & 3835 & *pawrhoij_unsym_all(at_indx)%rhoij_(indexkc+1:indexkc+cplex_eff,kspden) 3836 if (qphase==2) & 3837 & sumrho(1:cplex_eff,iafm,2)=sumrho(1:cplex_eff,iafm,2) & 3838 & +fact(1:cplex_eff)*zarot2 & 3839 & *pawrhoij_unsym_all(at_indx)%rhoij_(indexkc_q+1:indexkc_q+cplex_eff,kspden) 3840 end if 3841 3842 ! Rotate rhoij magnetization 3843 if (optrhoij==1.and.noncoll) then 3844 fact(1)=factsym(1)*factafm;fact(2)=factsym(2) 3845 do mu=1,3 3846 summag(1:cplex_eff,mu,1)=summag(1:cplex_eff,mu,1) & 3847 & +fact(1:cplex_eff)*zarot2 & 3848 & *pawrhoij_unsym_all(at_indx)%rhoij_(indexkc+1:indexkc+cplex_eff,1+mu) 3849 end do 3850 if (qphase==2) then 3851 do mu=1,3 3852 summag(1:cplex_eff,mu,2)=summag(1:cplex_eff,mu,2) & 3853 & +fact(1:cplex_eff)*zarot2 & 3854 & *pawrhoij_unsym_all(at_indx)%rhoij_(indexkc_q+1:indexkc_q+cplex_eff,1+mu) 3855 end do 3856 end if 3857 end if 3858 3859 ! Rotate gradients of rhoij 3860 if (choice>1) then 3861 fact(1)=factsym(1);fact(2)=factsym(2)*factafm !????? What? MT 3862 do mu=1,ngrhoij 3863 sumgr(1:cplex_eff,mu,1)=sumgr(1:cplex_eff,mu,1) & 3864 & +fact(1:cplex_eff)*zarot2*grad(mu,indexkc+1:indexkc+cplex_eff,kspden) 3865 end do 3866 if (qphase==2) then 3867 do mu=1,ngrhoij 3868 sumgr(1:cplex_eff,mu,2)=sumgr(1:cplex_eff,mu,2) & 3869 & +fact(1:cplex_eff)*zarot2*grad(mu,indexkc_q+1:indexkc_q+cplex_eff,kspden) 3870 end do 3871 end if 3872 end if 3873 3874 ! Rotate gradients of rhoij magnetization 3875 if (choice>1.and.noncoll) then 3876 fact(1)=factsym(1)*factafm;fact(2)=factsym(2) 3877 do mu=1,3 3878 do nu=1,ngrhoij 3879 summaggr(1:cplex_eff,nu,mu,1)=summaggr(1:cplex_eff,nu,mu,1) & 3880 & +fact(1:cplex_eff)*zarot2*grad(nu,indexkc+1:indexkc+cplex_eff,1+mu) 3881 end do 3882 end do 3883 if (qphase==2) then 3884 do mu=1,3 3885 do nu=1,ngrhoij 3886 summaggr(1:cplex_eff,nu,mu,2)=summaggr(1:cplex_eff,nu,mu,2) & 3887 & +fact(1:cplex_eff)*zarot2*grad(nu,indexkc_q+1:indexkc_q+cplex_eff,1+mu) 3888 end do 3889 end do 3890 end if 3891 end if 3892 3893 end do ! mi 3894 end do ! mj 3895 3896 ! Apply phase for phonons 3897 if (has_qphase) then 3898 !Remember, RHOij is stored as follows: 3899 ! RHOij= [rhoij(2klmn-1)+i.rhoij(2klmn)] 3900 ! +i.[rhoij(2lnm2_size+2klmn-1)+i.rhoij(2lmn2_size+2klmn)] 3901 if (optrhoij==1) then 3902 do iplex=1,cplex_rhoij 3903 rhoijc(1)=sumrho(iplex,iafm,1) 3904 rhoijc(2)=sumrho(iplex,iafm,2) 3905 sumrho(iplex,iafm,1)=phase(1)*rhoijc(1)-phase(2)*rhoijc(2) 3906 sumrho(iplex,iafm,2)=phase(1)*rhoijc(2)+phase(2)*rhoijc(1) 3907 end do 3908 if (noncoll) then 3909 do iplex=1,cplex_rhoij 3910 do mu=1,3 3911 rhoijc(1)=summag(iplex,mu,1) 3912 rhoijc(2)=summag(iplex,mu,2) 3913 summag(iplex,mu,1)=phase(1)*rhoijc(1)-phase(2)*rhoijc(2) 3914 summag(iplex,mu,2)=phase(1)*rhoijc(2)+phase(2)*rhoijc(1) 3915 end do 3916 end do 3917 end if 3918 end if 3919 if (choice>1) then 3920 do iplex=1,cplex_rhoij 3921 do mu=1,3 3922 rhoijc(1)=sumgr(iplex,mu,1) 3923 rhoijc(2)=sumgr(iplex,mu,2) 3924 sumgr(iplex,mu,1)=phase(1)*rhoijc(1)-phase(2)*rhoijc(2) 3925 sumgr(iplex,mu,2)=phase(1)*rhoijc(2)+phase(2)*rhoijc(1) 3926 end do 3927 end do 3928 if (noncoll) then 3929 do iplex=1,cplex_rhoij 3930 do mu=1,3 3931 do nu=1,ngrhoij 3932 rhoijc(1)=summaggr(iplex,nu,mu,1) 3933 rhoijc(2)=summaggr(iplex,nu,mu,2) 3934 summaggr(iplex,nu,mu,1)=phase(1)*rhoijc(1)-phase(2)*rhoijc(2) 3935 summaggr(iplex,nu,mu,2)=phase(1)*rhoijc(2)+phase(2)*rhoijc(1) 3936 end do 3937 end do 3938 end do 3939 end if 3940 end if 3941 end if 3942 3943 ! Add contribution of this rotation 3944 if (optrhoij==1) then 3945 do iq=1,qphase 3946 rotrho(1:cplex_eff,iafm,iq)=rotrho(1:cplex_eff,iafm,iq) & 3947 & +sumrho(1:cplex_eff,iafm,iq) 3948 end do 3949 end if 3950 3951 3952 ! Rotate vector fields in real space (forces, magnetization, etc...) 3953 ! Should use symrel^-1 but use transpose[symrec] instead 3954 ! ===== Rhoij magnetization ==== 3955 if (noncoll.and.optrhoij==1) then 3956 do iq=1,qphase 3957 do nu=1,3 3958 do mu=1,3 3959 rotmag(1:cplex_eff,mu,iq)=rotmag(1:cplex_eff,mu,iq) & 3960 & +symrec_cart(mu,nu,irot)*summag(1:cplex_eff,nu,iq) 3961 end do 3962 end do 3963 end do 3964 end if 3965 ! ===== Derivatives vs atomic positions ==== 3966 if (choice==2.or.choice==23.or.choice==24) then 3967 do iq=1,qphase 3968 do nu=1,3 3969 nushift=nu+ishift2 3970 do mu=1,3 3971 mushift=mu+ishift2 3972 rotgr(1:cplex_eff,mushift,iafm,iq)=rotgr(1:cplex_eff,mushift,iafm,iq) & 3973 & +dble(symrec(mu,nu,irot))*sumgr(1:cplex_eff,nushift,iq) 3974 end do 3975 end do 3976 end do 3977 if (noncoll) then 3978 do iq=1,qphase 3979 do mub=1,3 ! Loop on magnetization components 3980 do mua=1,3 ! Loop on gradients 3981 mushift=mua+ishift2 3982 sum1(:)=zero;xsym(1:3)=dble(symrec(mua,1:3,irot)) 3983 do nu=1,3 3984 syma=symrec_cart(mub,nu,irot) 3985 sum1(1:cplex_eff)=sum1(1:cplex_eff)+syma & 3986 & *(summaggr(1:cplex_eff,ishift2+1,nu,iq)*xsym(1) & 3987 & +summaggr(1:cplex_eff,ishift2+2,nu,iq)*xsym(2) & 3988 & +summaggr(1:cplex_eff,ishift2+3,nu,iq)*xsym(3)) 3989 end do 3990 rotmaggr(1:cplex_eff,mushift,mub,iq)= & 3991 & rotmaggr(1:cplex_eff,mushift,mub,iq)+sum1(1:cplex_eff) 3992 end do 3993 end do 3994 end do 3995 end if 3996 end if 3997 ! ===== Derivatives vs strain ==== 3998 if (choice==3.or.choice==23) then 3999 do iq=1,qphase 4000 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) 4001 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) 4002 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) 4003 work1(1:cplex_eff,3,1)=work1(1:cplex_eff,1,3);work1(1:cplex_eff,3,2)=work1(1:cplex_eff,2,3) 4004 work1(1:cplex_eff,2,1)=work1(1:cplex_eff,1,2) 4005 do mu=1,6 4006 mushift=mu+ishift3 4007 mua=alpha(mu);mub=beta(mu) 4008 sum1(:)=zero;xsym(1:3)=dble(symrec(mub,1:3,irot)) 4009 do nu=1,3 4010 syma=dble(symrec(mua,nu,irot)) 4011 sum1(1:cplex_eff)=sum1(1:cplex_eff) & 4012 & +syma*(work1(1:cplex_eff,nu,1)*xsym(1) & 4013 & +work1(1:cplex_eff,nu,2)*xsym(2) & 4014 & +work1(1:cplex_eff,nu,3)*xsym(3)) 4015 end do 4016 rotgr(1:cplex_eff,mushift,iafm,iq)= & 4017 & rotgr(1:cplex_eff,mushift,iafm,iq)+sum1(1:cplex_eff) 4018 end do 4019 end do 4020 end if 4021 ! ===== Second derivatives vs atomic positions ==== 4022 if (choice==4.or.choice==24) then 4023 do iq=1,qphase 4024 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) 4025 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) 4026 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) 4027 work1(1:cplex_eff,3,1)=work1(1:cplex_eff,1,3);work1(1:cplex_eff,3,2)=work1(1:cplex_eff,2,3) 4028 work1(1:cplex_eff,2,1)=work1(1:cplex_eff,1,2) 4029 do mu=1,6 4030 mushift=mu+ishift4 4031 mua=alpha(mu);mub=beta(mu) 4032 sum1(:)=zero 4033 xsym(1:3)=dble(symrec(mub,1:3,irot)) 4034 do nu=1,3 4035 syma=dble(symrec(mua,nu,irot)) 4036 sum1(1:cplex_eff)=sum1(1:cplex_eff) & 4037 & +syma*(work1(1:cplex_eff,nu,1)*xsym(1) & 4038 & +work1(1:cplex_eff,nu,2)*xsym(2) & 4039 & +work1(1:cplex_eff,nu,3)*xsym(3)) 4040 end do 4041 rotgr(1:cplex_eff,mushift,iafm,iq)= & 4042 & rotgr(1:cplex_eff,mushift,iafm,iq)+sum1(1:cplex_eff) 4043 end do 4044 end do 4045 end if 4046 4047 end do ! End loop over symmetries 4048 4049 4050 ! Store average result (over symmetries) 4051 ! -------------------------------------- 4052 4053 ! Rhoij 4054 if (optrhoij==1) then 4055 do iq=1,qphase 4056 klmn1q=klmn1+(iq-1)*lmn2_size*cplex_rhoij 4057 pawrhoij(iatm)%rhoijp(klmn1q,ispden)=rotrho(1,1,iq)/nsym_used(1) 4058 if (cplex_rhoij==2) then 4059 if (cplex_eff==1) ro=pawrhoij_unsym_all(iatom)%rhoij_(klmn1q+1,ispden) 4060 if (cplex_eff==2) ro=rotrho(2,1,iq)/nsym_used(1) 4061 pawrhoij(iatm)%rhoijp(klmn1q+1,ispden)=ro 4062 end if 4063 end do 4064 end if 4065 4066 ! Rhoij magnetization 4067 if (noncoll.and.optrhoij==1) then 4068 do mu=2,4 4069 do iq=1,qphase 4070 klmn1q=klmn1+(iq-1)*lmn2_size*cplex_rhoij 4071 if (use_zeromag_) then 4072 pawrhoij(iatm)%rhoijp(klmn1q,mu)=zero 4073 else 4074 pawrhoij(iatm)%rhoijp(klmn1q,mu)=rotmag(1,mu-1,iq)/nsym_used(1) 4075 end if 4076 if (cplex_rhoij==2) then 4077 if (cplex_eff==1) ro=pawrhoij_unsym_all(iatom)%rhoij_(klmn1q+1,mu) 4078 if (cplex_eff==2) ro=rotmag(2,mu-1,iq)/nsym_used(1) 4079 pawrhoij(iatm)%rhoijp(klmn1q+1,mu)=ro 4080 end if 4081 end do 4082 end do 4083 end if 4084 4085 ! Rhoij^down when antiferro 4086 if (antiferro.and.optrhoij==1) then 4087 if (nsym_used(2)>0) then 4088 do iq=1,qphase 4089 klmn1q=klmn1+(iq-1)*lmn2_size*cplex_rhoij 4090 pawrhoij(iatm)%rhoijp(klmn1q,2)=rotrho(1,2,iq)/nsym_used(2) 4091 if (cplex_rhoij==2) then 4092 if (cplex_eff==1) ro=pawrhoij_unsym_all(iatom)%rhoij_(klmn1q+1,2) 4093 if (cplex_eff==2) ro=rotrho(2,2,iq)/nsym_used(2) 4094 pawrhoij(iatm)%rhoijp(klmn1q+1,2)=ro 4095 end if 4096 end do 4097 end if 4098 end if 4099 4100 ! Gradients of rhoij 4101 if (choice>1) then 4102 do iq=1,qphase 4103 klmn1q=klmn1+(iq-1)*lmn2_size*cplex_rhoij 4104 do iplex=1,cplex_eff 4105 do mu=1,ngrhoij 4106 pawrhoij(iatm)%grhoij(mu,klmn1q,ispden)=rotgr(iplex,mu,1,iq)/nsym_used(1) 4107 end do 4108 if (noncoll) then 4109 if (use_zeromag_.and.iplex==1) then 4110 pawrhoij(iatm)%grhoij(mu,klmn1q,2:4)=zero 4111 else 4112 do nu=1,3 4113 pawrhoij(iatm)%grhoij(mu,klmn1q,1+nu)=rotmaggr(iplex,mu,nu,iq)/nsym_used(1) 4114 end do 4115 end if 4116 end if 4117 if (antiferro.and.nsym_used(2)>0) then 4118 do mu=1,ngrhoij 4119 pawrhoij(iatm)%grhoij(mu,klmn1q,ispden)=rotgr(iplex,mu,2,iq)/nsym_used(2) 4120 end do 4121 end if 4122 klmn1q=klmn1q+1 4123 end do 4124 !if cplex_eff<cplex_rhoij, imaginary part of grhoij is unchanged 4125 end do 4126 end if 4127 4128 4129 il0=il;iln0=iln ! End loops over (il,im) and (jl,jm) 4130 end do 4131 jl0=jl;jln0=jln 4132 end do 4133 4134 end do ! End loop over ispden 4135 4136 ! Select non-zero elements of rhoij 4137 if (optrhoij==1) then 4138 call pawrhoij_filter(pawrhoij(iatm)%rhoijp,pawrhoij(iatm)%rhoijselect,& 4139 & pawrhoij(iatm)%nrhoijsel,cplex_rhoij,qphase,lmn2_size,& 4140 & pawrhoij(iatm)%nspden) 4141 end if 4142 4143 ! Add new rhoij to rhoij residual 4144 if (optrhoij==1.and.use_res) then 4145 do ispden=1,pawrhoij(iatm)%nspden 4146 do iq=1,qphase 4147 iq0=(iq-1)*lmn2_size*cplex_rhoij 4148 if (cplex_rhoij==1) then 4149 do irhoij=1,pawrhoij(iatm)%nrhoijsel 4150 klmn1=iq0+pawrhoij(iatm)%rhoijselect(irhoij) ; jrhoij=iq0+irhoij 4151 pawrhoij(iatm)%rhoijres(klmn1,ispden)= & 4152 & pawrhoij(iatm)%rhoijres(klmn1,ispden) & 4153 & +pawrhoij(iatm)%rhoijp(jrhoij,ispden) 4154 end do 4155 else 4156 do irhoij=1,pawrhoij(iatm)%nrhoijsel 4157 klmn1=iq0+2*pawrhoij(iatm)%rhoijselect(irhoij) ; jrhoij=iq0+2*irhoij 4158 pawrhoij(iatm)%rhoijres(klmn1-1:klmn1,ispden)= & 4159 & pawrhoij(iatm)%rhoijres(klmn1-1:klmn1,ispden) & 4160 & +pawrhoij(iatm)%rhoijp(jrhoij-1:jrhoij,ispden) 4161 end do 4162 end if 4163 end do 4164 end do 4165 end if 4166 4167 end do ! End loop over atoms 4168 4169 if (noncoll) then 4170 LIBPAW_DEALLOCATE(symrec_cart) 4171 end if 4172 if (choice>1) then 4173 if (.not.paral_atom_unsym) then 4174 do iatm=1,nrhoij 4175 LIBPAW_DEALLOCATE(tmp_grhoij(iatm)%value) 4176 end do 4177 LIBPAW_DATATYPE_DEALLOCATE(tmp_grhoij) 4178 end if 4179 LIBPAW_DEALLOCATE(sumgr) 4180 LIBPAW_DEALLOCATE(rotgr) 4181 if (noncoll) then 4182 LIBPAW_DEALLOCATE(summaggr) 4183 LIBPAW_DEALLOCATE(rotmaggr) 4184 end if 4185 end if 4186 if(paral_atom_unsym) then 4187 call pawrhoij_free(pawrhoij_unsym_all) 4188 LIBPAW_DATATYPE_DEALLOCATE(pawrhoij_unsym_all) 4189 end if 4190 4191 4192 else ! nsym>1 4193 4194 ! ********************************************************************* 4195 ! If nsym==1, only copy rhoij_ into rhoij 4196 ! also has to fill rhoijselect array 4197 4198 if (antiferro) then 4199 msg=' In the antiferromagnetic case, nsym cannot be 1' 4200 LIBPAW_BUG(msg) 4201 end if 4202 4203 if (optrhoij==1) then 4204 4205 do iatm=1,nrhoij 4206 iatom=iatm;if ((paral_atom).and.(.not.paral_atom_unsym)) iatom=my_atmtab(iatm) 4207 cplex_rhoij=pawrhoij(iatm)%cplex_rhoij 4208 qphase=pawrhoij(iatm)%qphase 4209 lmn2_size=pawrhoij(iatm)%lmn2_size 4210 use_res=(pawrhoij(iatm)%use_rhoijres>0) 4211 4212 ! Store -rhoij_input in rhoij residual 4213 if (use_res) then 4214 pawrhoij(iatm)%rhoijres(:,:)=zero 4215 do iq=1,qphase 4216 iq0=(iq-1)*lmn2_size*cplex_rhoij 4217 if (cplex_rhoij==1) then 4218 do ispden=1,pawrhoij(iatm)%nspden 4219 do irhoij=1,pawrhoij(iatm)%nrhoijsel 4220 klmn=iq0+pawrhoij(iatm)%rhoijselect(irhoij);jrhoij=iq0+irhoij 4221 pawrhoij(iatm)%rhoijres(klmn,ispden)=-pawrhoij(iatm)%rhoijp(jrhoij,ispden) 4222 end do 4223 end do 4224 else 4225 do ispden=1,pawrhoij(iatm)%nspden 4226 do irhoij=1,pawrhoij(iatm)%nrhoijsel 4227 klmn1=iq0+2*pawrhoij(iatm)%rhoijselect(irhoij);jrhoij=iq0+2*irhoij 4228 pawrhoij(iatm)%rhoijres(klmn1-1,ispden)=-pawrhoij(iatm)%rhoijp(jrhoij-1,ispden) 4229 pawrhoij(iatm)%rhoijres(klmn1 ,ispden)=-pawrhoij(iatm)%rhoijp(jrhoij ,ispden) 4230 end do 4231 end do 4232 end if 4233 end do 4234 end if 4235 4236 ! Select non-zero elements of input rhoij 4237 call pawrhoij_filter(pawrhoij(iatm)%rhoijp,pawrhoij(iatm)%rhoijselect,& 4238 & pawrhoij(iatm)%nrhoijsel,cplex_rhoij,qphase,lmn2_size,& 4239 & pawrhoij(iatm)%nspden,rhoij_input=pawrhoij_unsym(iatom)%rhoij_) 4240 4241 ! Add new rhoij to rhoij residual 4242 if (use_res) then 4243 do ispden=1,pawrhoij(iatm)%nspden 4244 do iq=1,qphase 4245 iq0=(iq-1)*lmn2_size*cplex_rhoij 4246 if (cplex_rhoij==1) then 4247 do irhoij=1,pawrhoij(iatm)%nrhoijsel 4248 klmn1=iq0+pawrhoij(iatm)%rhoijselect(irhoij) ; jrhoij=iq0+irhoij 4249 pawrhoij(iatm)%rhoijres(klmn1,ispden)= & 4250 & pawrhoij(iatm)%rhoijres(klmn1,ispden) & 4251 & +pawrhoij(iatm)%rhoijp(jrhoij,ispden) 4252 end do 4253 else 4254 do irhoij=1,pawrhoij(iatm)%nrhoijsel 4255 klmn1=iq0+2*pawrhoij(iatm)%rhoijselect(irhoij) ; jrhoij=iq0+2*irhoij 4256 pawrhoij(iatm)%rhoijres(klmn1-1:klmn1,ispden)= & 4257 & pawrhoij(iatm)%rhoijres(klmn1-1:klmn1,ispden) & 4258 & +pawrhoij(iatm)%rhoijp(jrhoij-1:jrhoij,ispden) 4259 end do 4260 end if 4261 end do 4262 end do 4263 end if 4264 4265 end do ! iatm 4266 end if ! optrhoij 4267 4268 end if 4269 4270 !********************************************************************* 4271 !Printing of symetrized Rhoij 4272 if (nrhoij>0.and.optrhoij==1.and.pawprtvol/=-10001) then 4273 wrt_mode='COLL';if (paral_atom) wrt_mode='PERS' 4274 pertstrg="RHOIJ";if (ipert>0) pertstrg="RHOIJ(1)" 4275 natinc=1;if(nrhoij>1.and.pawprtvol>=0) natinc=nrhoij-1 4276 write(msg, '(7a)') ch10," PAW TEST:",ch10,& 4277 & ' ========= Values of ',trim(pertstrg),' after symetrization =========',ch10 4278 call wrtout(std_out,msg,wrt_mode) 4279 do iatm=1,nrhoij,natinc 4280 iatom=iatm; if (paral_atom) iatom=my_atmtab(iatm) 4281 if (nrhoij==1.and.ipert>0.and.ipert<=natom) iatom=ipert 4282 call pawrhoij_print_rhoij(pawrhoij(iatm)%rhoijp,pawrhoij(iatm)%cplex_rhoij,& 4283 & pawrhoij(iatm)%qphase,iatom,natom,& 4284 & rhoijselect=pawrhoij(iatm)%rhoijselect,unit=std_out,& 4285 & opt_prtvol=pawprtvol,mode_paral=wrt_mode) 4286 end do 4287 call wrtout(std_out,"",wrt_mode) 4288 end if 4289 4290 !Destroy atom table used for parallelism 4291 call free_my_atmtab(my_atmtab,my_atmtab_allocated) 4292 4293 !********************************************************************* 4294 !Small function: convert a symmetry operation 4295 !from reduced coordinates (integers) to cartesian coordinates (reals) 4296 contains 4297 function symrhoij_symcart(aprim,bprim,symred) 4298 4299 real(dp) :: symrhoij_symcart(3,3) 4300 integer,intent(in) :: symred(3,3) 4301 real(dp),intent(in) :: aprim(3,3),bprim(3,3) 4302 integer :: ii,jj,kk 4303 real(dp) :: tmp(3,3) 4304 symrhoij_symcart=zero;tmp=zero 4305 do kk=1,3 4306 do jj=1,3 4307 do ii=1,3 4308 tmp(ii,jj)=tmp(ii,jj)+bprim(ii,kk)*dble(symred(jj,kk)) 4309 end do 4310 end do 4311 end do 4312 do kk=1,3 4313 do jj=1,3 4314 do ii=1,3 4315 symrhoij_symcart(ii,jj)=symrhoij_symcart(ii,jj)+aprim(ii,kk)*tmp(jj,kk) 4316 end do 4317 end do 4318 end do 4319 end function symrhoij_symcart 4320 4321 end subroutine pawrhoij_symrhoij
m_pawrhoij/pawrhoij_type [ Types ]
[ Top ] [ m_pawrhoij ] [ Types ]
NAME
pawrhoij_type
FUNCTION
This structured datatype contains rhoij quantities (occucpancies) and related data, used in PAW calculations.
SOURCE
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
m_pawrhoij/pawrhoij_unpack [ Functions ]
[ Top ] [ m_pawrhoij ] [ Functions ]
NAME
pawrhoij_unpack
FUNCTION
Unpack the values store in rhoijp copying them to the rhoij_ array.
SIDE EFFECTS
rhoij(:) <pawrhoij_type)>= input/output datastructure * In output the rhoij_ array is filled with the values stored in the packed array rhoijp. * If use_rhoij_/=1, rhoij_ is allocated and the corresponding flag is set to 1.
SOURCE
2734 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 natom = SIZE(rhoij) ; if (natom==0) return 2748 2749 do iat=1,natom 2750 2751 lmn2_size =rhoij(iat)%lmn2_size 2752 cplex_rhoij = rhoij(iat)%cplex_rhoij 2753 qphase = rhoij(iat)%qphase 2754 nsp2 = rhoij(iat)%nsppol;if (rhoij(iat)%nspden==4) nsp2=4 2755 2756 if (rhoij(iat)%use_rhoij_/=1) then ! Have to allocate rhoij 2757 LIBPAW_ALLOCATE(rhoij(iat)%rhoij_,(cplex_rhoij*qphase*lmn2_size,nsp2)) 2758 rhoij(iat)%use_rhoij_=1 2759 end if 2760 rhoij(iat)%rhoij_ = zero 2761 2762 do isppol=1,nsp2 2763 do iphase=1,qphase 2764 i0=(iphase-1)*lmn2_size*cplex_rhoij 2765 if (cplex_rhoij==1) then 2766 do isel=1,rhoij(iat)%nrhoijsel ! Looping over non-zero ij elements. 2767 klmn = rhoij(iat)%rhoijselect(isel) 2768 rhoij(iat)%rhoij_(i0+klmn,isppol) = rhoij(iat)%rhoijp(i0+isel,isppol) 2769 end do 2770 else 2771 do isel=1,rhoij(iat)%nrhoijsel ! Looping over non-zero ij elements. 2772 klmn = rhoij(iat)%rhoijselect(isel) 2773 rhoij(iat)%rhoij_(i0+2*klmn-1:i0+2*klmn,isppol) = & 2774 & 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