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
3269 subroutine pawrhoij_print_rhoij(rhoij,cplex,qphase,iatom,natom,& 3270 & rhoijselect,test_value,title_msg,unit,opt_prtvol,& 3271 & l_only,indlmn,mode_paral) ! Optional arguments 3272 3273 !Arguments ------------------------------------ 3274 !scalars 3275 integer,intent(in) :: cplex,qphase,iatom,natom 3276 integer,optional,intent(in) :: opt_prtvol,l_only,unit 3277 real(dp),intent(in),optional :: test_value 3278 character(len=4),optional,intent(in) :: mode_paral 3279 character(len=100),optional,intent(in) :: title_msg 3280 !arrays 3281 integer,optional,intent(in) :: indlmn(:,:) 3282 integer,optional,intent(in),target :: rhoijselect(:) 3283 real(dp),intent(in),target :: rhoij(:,:) 3284 3285 !Local variables------------------------------- 3286 !scalars 3287 character(len=8),parameter :: dspin(6)=(/"up ","down ","dens (n)","magn (x)","magn (y)","magn (z)"/) 3288 integer :: irhoij,kk,my_cplex,my_lmn_size,my_lmn2_size,my_l_only,my_nspden 3289 integer :: my_opt_pack,my_opt_sym,my_prtvol,my_unt,nrhoijsel,rhoij_size 3290 real(dp) :: my_test_value,test_value_eff 3291 character(len=4) :: my_mode 3292 character(len=2000) :: msg 3293 !arrays 3294 integer,target :: idum(0) 3295 integer,pointer :: l_index(:),my_rhoijselect(:) 3296 real(dp),pointer :: rhoij_(:) 3297 3298 ! ************************************************************************* 3299 3300 !Optional arguments 3301 my_unt =std_out ; if (PRESENT(unit )) my_unt =unit 3302 my_mode ='COLL' ; if (PRESENT(mode_paral)) my_mode =mode_paral 3303 my_prtvol=1 ; if (PRESENT(opt_prtvol)) my_prtvol=opt_prtvol 3304 my_test_value=-one; if (PRESENT(test_value)) my_test_value=test_value 3305 my_l_only=-1 ; if (PRESENT(l_only)) my_l_only=l_only 3306 3307 if (my_l_only>=0.and.(.not.present(indlmn))) then 3308 msg='pawrhoij_print_rhoij: l_only>=0 and indlmn not present!' 3309 LIBPAW_BUG(msg) 3310 end if 3311 3312 !Title 3313 if (present(title_msg)) then 3314 if (trim(title_msg)/='') then 3315 write(msg, '(2a)') ch10,trim(title_msg) 3316 call wrtout(my_unt,msg,my_mode) 3317 end if 3318 end if 3319 3320 !Inits 3321 my_nspden=size(rhoij,2) 3322 my_lmn2_size=size(rhoij,1)/cplex/qphase 3323 my_lmn_size=int(dsqrt(two*dble(my_lmn2_size))) 3324 my_cplex=merge(cplex,2,qphase==1) 3325 my_opt_sym=2 3326 3327 !Packed storage 3328 my_opt_pack=0 3329 rhoij_size=my_lmn2_size 3330 my_rhoijselect => idum 3331 if (PRESENT(rhoijselect)) then 3332 nrhoijsel=count(rhoijselect(:)>0) 3333 if (nrhoijsel>0) then 3334 my_opt_pack=1 3335 rhoij_size=nrhoijsel 3336 my_rhoijselect => rhoijselect(1:nrhoijsel) 3337 end if 3338 end if 3339 3340 if (my_l_only<0) then 3341 l_index => idum 3342 else 3343 LIBPAW_POINTER_ALLOCATE(l_index,(my_lmn_size)) 3344 do kk=1,my_lmn_size 3345 l_index(kk)=indlmn(1,kk) 3346 end do 3347 end if 3348 3349 if (qphase==2) then 3350 LIBPAW_DATATYPE_ALLOCATE(rhoij_,(2*rhoij_size)) 3351 end if 3352 3353 ! === Loop over Rho_ij components === 3354 do irhoij=1,my_nspden 3355 3356 !Rebuild rhoij according to qphase 3357 if (qphase==1) then 3358 rhoij_ => rhoij(1:cplex*rhoij_size,irhoij) 3359 else 3360 if (cplex==1) then 3361 do kk=1,rhoij_size 3362 rhoij_(2*kk-1)=rhoij(kk,irhoij) 3363 rhoij_(2*kk )=rhoij(kk+my_lmn2_size,irhoij) 3364 end do 3365 else 3366 do kk=1,rhoij_size 3367 rhoij_(2*kk-1)=rhoij(2*kk-1,irhoij)-rhoij(2*kk +2*my_lmn2_size,irhoij) 3368 rhoij_(2*kk )=rhoij(2*kk ,irhoij)+rhoij(2*kk-1+2*my_lmn2_size,irhoij) 3369 end do 3370 end if 3371 end if 3372 3373 !Subtitle 3374 if (natom>1.or.my_nspden>1) then 3375 if (my_l_only<0) then 3376 if (my_nspden==1) write(msg,'(a,i3)') ' Atom #',iatom 3377 if (my_nspden==2) write(msg,'(a,i3,a,i1)')' Atom #',iatom,' - Spin component ',irhoij 3378 if (my_nspden==4) write(msg,'(a,i3,2a)') ' Atom #',iatom,' - Component ',trim(dspin(irhoij+2*(my_nspden/4))) 3379 else 3380 if (my_nspden==1) write(msg,'(a,i3,a,i1,a)') ' Atom #',iatom,& 3381 & ' - L=',my_l_only,' ONLY' 3382 if (my_nspden==2) write(msg,'(a,i3,a,i1,a,i1)')' Atom #',iatom,& 3383 & ' - L=',my_l_only,' ONLY - Spin component ',irhoij 3384 if (my_nspden==4) write(msg,'(a,i3,a,i1,3a)') ' Atom #',iatom,& 3385 & ' - L=',my_l_only,' ONLY - Component ',trim(dspin(irhoij+2*(my_nspden/4))) 3386 end if 3387 call wrtout(my_unt,msg,my_mode) 3388 else if (my_l_only>=0) then 3389 write(msg,'(a,i1,a)') ' L=',my_l_only,' ONLY' 3390 call wrtout(my_unt,msg,my_mode) 3391 end if 3392 3393 !Printing 3394 test_value_eff=-one;if(my_test_value>zero.and.irhoij==1) test_value_eff=my_test_value 3395 call pawio_print_ij(my_unt,rhoij_,rhoij_size,my_cplex,my_lmn_size,my_l_only,l_index,my_opt_pack,& 3396 & my_prtvol,my_rhoijselect,test_value_eff,1,opt_sym=my_opt_sym,& 3397 & mode_paral=my_mode) 3398 3399 end do !irhoij 3400 3401 if (qphase==2) then 3402 LIBPAW_DATATYPE_DEALLOCATE(rhoij_) 3403 end if 3404 if (my_l_only>=0) then 3405 LIBPAW_POINTER_DEALLOCATE(l_index) 3406 end if 3407 3408 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
3067 subroutine pawrhoij_filter(rhoij,rhoijselect,nselect,cplex,qphase,lmn2_size,nspden, & 3068 & rhoij_input) ! optional argument 3069 3070 !Arguments ------------------------------------ 3071 !scalars 3072 integer,intent(in) :: lmn2_size,cplex,qphase,nspden 3073 integer,intent(out) :: nselect 3074 !arrays 3075 integer,intent(out) :: rhoijselect(lmn2_size) 3076 real(dp),intent(inout),target :: rhoij(cplex*qphase*lmn2_size,nspden) 3077 real(dp),intent(in),optional,target :: rhoij_input(cplex*qphase*lmn2_size,nspden) 3078 3079 !Local variables------------------------------- 3080 !scalars 3081 real(dp),parameter :: tol_rhoij=tol10 3082 integer :: isp,klmn,klmn1,klmn2,nsel1,nsel2 3083 !arrays 3084 real(dp),pointer :: rhoij_in(:,:) 3085 3086 ! ************************************************************************* 3087 3088 nselect=0 3089 rhoijselect(:)=0 3090 3091 if (present(rhoij_input)) then 3092 rhoij_in => rhoij_input 3093 else 3094 rhoij_in => rhoij 3095 end if 3096 3097 !Treat each cplex/qphase case separately 3098 3099 if (cplex==1) then 3100 3101 if (qphase==1) then 3102 3103 do klmn=1,lmn2_size 3104 if (any(abs(rhoij_in(klmn,:))>tol_rhoij)) then 3105 nselect=nselect+1 3106 rhoijselect(nselect)=klmn 3107 do isp=1,nspden 3108 rhoij(nselect,isp)=rhoij_in(klmn,isp) 3109 end do 3110 end if 3111 end do 3112 3113 else if (qphase==2) then 3114 3115 do klmn=1,lmn2_size 3116 klmn2=klmn+lmn2_size 3117 if (any(abs(rhoij_in(klmn,:))>tol_rhoij).or. & 3118 & any(abs(rhoij_in(klmn2,:))>tol_rhoij)) then 3119 nselect=nselect+1 ; nsel2=nselect+lmn2_size 3120 rhoijselect(nselect)=klmn 3121 do isp=1,nspden 3122 rhoij(nselect,isp)=rhoij_in(klmn ,isp) 3123 rhoij(nsel2 ,isp)=rhoij_in(klmn2,isp) 3124 end do 3125 end if 3126 end do 3127 3128 end if 3129 3130 else ! cplex=2 3131 3132 if (qphase==1) then 3133 do klmn=1,lmn2_size 3134 klmn1=2*klmn 3135 if (any(abs(rhoij_in(klmn1-1:klmn1,:))>tol_rhoij)) then 3136 nselect=nselect+1 ; nsel1=2*nselect 3137 rhoijselect(nselect)=klmn 3138 do isp=1,nspden 3139 rhoij(nsel1-1,isp)=rhoij_in(klmn1-1,isp) 3140 rhoij(nsel1 ,isp)=rhoij_in(klmn1 ,isp) 3141 end do 3142 end if 3143 end do 3144 3145 else if (qphase==2) then 3146 3147 do klmn=1,lmn2_size 3148 klmn1=2*klmn ; klmn2=klmn1+lmn2_size 3149 if (any(abs(rhoij_in(klmn1-1:klmn1,:))>tol_rhoij).or. & 3150 & any(abs(rhoij_in(klmn2-1:klmn2,:))>tol_rhoij)) then 3151 nselect=nselect+1 ; nsel1=2*nselect ; nsel2=nsel1+lmn2_size 3152 rhoijselect(nselect)=klmn 3153 do isp=1,nspden 3154 rhoij(nsel1-1,isp)=rhoij_in(klmn1-1,isp) 3155 rhoij(nsel1 ,isp)=rhoij_in(klmn1 ,isp) 3156 rhoij(nsel2-1,isp)=rhoij_in(klmn2-1,isp) 3157 rhoij(nsel2 ,isp)=rhoij_in(klmn2 ,isp) 3158 end do 3159 end if 3160 end do 3161 3162 end if 3163 endif 3164 3165 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
2832 subroutine pawrhoij_free_unpacked(rhoij) 2833 2834 !Arguments ------------------------------------ 2835 !scalars 2836 !arrays 2837 type(pawrhoij_type),intent(inout) :: rhoij(:) 2838 2839 !Local variables------------------------------- 2840 integer :: iat,nrhoij 2841 2842 ! ************************************************************************* 2843 2844 nrhoij=SIZE(rhoij);if (nrhoij==0) return 2845 2846 do iat=1,nrhoij 2847 2848 if (allocated(rhoij(iat)%rhoij_)) then 2849 LIBPAW_DEALLOCATE(rhoij(iat)%rhoij_) 2850 end if 2851 rhoij(iat)%use_rhoij_=0 2852 2853 end do 2854 2855 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
2784 subroutine pawrhoij_init_unpacked(rhoij) 2785 2786 !Arguments ------------------------------------ 2787 !scalars 2788 !arrays 2789 type(pawrhoij_type),intent(inout) :: rhoij(:) 2790 2791 !Local variables------------------------------- 2792 integer :: cplex_rhoij,iat,lmn2_size,nrhoij,nsp2,qphase 2793 2794 ! ************************************************************************* 2795 2796 nrhoij=SIZE(rhoij);if (nrhoij==0) return 2797 2798 do iat=1,nrhoij 2799 2800 lmn2_size =rhoij(iat)%lmn2_size 2801 cplex_rhoij = rhoij(iat)%cplex_rhoij 2802 qphase = rhoij(iat)%qphase 2803 nsp2 = rhoij(iat)%nsppol;if (rhoij(iat)%nspden==4) nsp2=4 2804 2805 if (allocated(rhoij(iat)%rhoij_)) then 2806 LIBPAW_DEALLOCATE(rhoij(iat)%rhoij_) 2807 end if 2808 LIBPAW_ALLOCATE(rhoij(iat)%rhoij_,(cplex_rhoij*qphase*lmn2_size,nsp2)) 2809 rhoij(iat)%use_rhoij_=1 2810 rhoij(iat)%rhoij_=zero 2811 2812 end do 2813 2814 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
3192 subroutine pawrhoij_inquire_dim(cplex,cpxocc,nspden,qpt,spnorb, & 3193 & cplex_rhoij,qphase_rhoij,nspden_rhoij) 3194 3195 !Arguments --------------------------------------------- 3196 !scalars 3197 integer,optional,intent(in) :: cplex,cpxocc,nspden,spnorb 3198 integer,optional,intent(out) :: cplex_rhoij,qphase_rhoij,nspden_rhoij 3199 !arrays 3200 real(dp),optional,intent(in) :: qpt(3) 3201 3202 !Local variables --------------------------------------- 3203 character(len=100) :: msg 3204 3205 !************************************************************************ 3206 3207 !cplex_rhoij 3208 if (present(cplex_rhoij)) then 3209 cplex_rhoij=1 3210 if (present(cpxocc)) cplex_rhoij=max(cplex_rhoij,cpxocc) 3211 end if 3212 3213 !qphase_rhoij 3214 if (present(qphase_rhoij)) then 3215 qphase_rhoij=1 3216 if (present(cplex).and.present(qpt)) then 3217 msg='only one argument cplex or qpt should be passed!' 3218 LIBPAW_BUG(msg) 3219 end if 3220 if (present(cplex)) qphase_rhoij=merge(1,2,cplex==1) 3221 if (present(qpt)) then 3222 if (any(abs(qpt(:))>tol8)) qphase_rhoij=2 3223 end if 3224 end if 3225 3226 !nspden_rhoij 3227 if (present(nspden_rhoij)) then 3228 nspden_rhoij=1 3229 if (present(nspden)) nspden_rhoij=nspden 3230 if (present(spnorb)) nspden_rhoij=merge(nspden_rhoij,4,spnorb<=0) 3231 end if 3232 3233 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 logical :: paral_atom 2302 character(len=500) :: msg 2303 !arrays 2304 integer,allocatable :: ibuffer(:),nsel44(:,:),nsel56(:) 2305 integer,pointer :: my_atmtab(:) 2306 real(dp), allocatable :: buffer(:) 2307 real(dp),pointer :: rhoij_tmp(:) 2308 2309 ! ************************************************************************* 2310 2311 my_natom=SIZE(pawrhoij);if (my_natom==0) return 2312 my_nspden=nspden_in 2313 natom=size(typat) 2314 paral_atom=(my_natom/=natom) 2315 if (present(mpi_atmtab)) then 2316 if (.not.associated(mpi_atmtab)) then 2317 msg='mpi_atmtab not associated (pawrhoij_io)' 2318 LIBPAW_BUG(msg) 2319 end if 2320 my_atmtab=>mpi_atmtab 2321 else if (my_natom/=natom) then 2322 msg='my_natom /=natom, mpi_atmtab should be in argument (pawrhoij_io)' 2323 LIBPAW_BUG(msg) 2324 end if 2325 2326 iomode = fort_binary 2327 if (PRESENT(form)) then 2328 select case (libpaw_to_upper(form)) 2329 case ("FORMATTED") 2330 iomode = fort_formatted 2331 case ("NETCDF") 2332 iomode = netcdf_io 2333 case default 2334 LIBPAW_ERROR("Wrong form: "//trim(form)) 2335 end select 2336 end if 2337 2338 #ifndef LIBPAW_HAVE_NETCDF 2339 if (iomode == netcdf_io) then 2340 LIBPAW_ERROR("iomode == netcdf_io but netcdf library is missing.") 2341 end if 2342 #endif 2343 ncid = unitfi 2344 2345 select case (rdwr_mode(1:1)) 2346 2347 case ("R","r") ! Reading the Rhoij tab 2348 2349 if ((headform>=44).and.(headform<56)) then 2350 LIBPAW_ALLOCATE(nsel44,(nspden_in,natom)) 2351 if (iomode == fort_binary) then 2352 read(unitfi ) ((nsel44(ispden,iatom),ispden=1,nspden_in),iatom=1,natom) 2353 else if (iomode == fort_formatted) then 2354 read(unitfi,*) ((nsel44(ispden,iatom),ispden=1,nspden_in),iatom=1,natom) 2355 #ifdef LIBPAW_HAVE_NETCDF 2356 else if (iomode == netcdf_io) then 2357 LIBPAW_ERROR("header in 44-56 not compatible with Netcdf") 2358 #endif 2359 end if 2360 call pawrhoij_alloc(pawrhoij,1,nspden_in,nspinor_in,nsppol_in,typat,lmnsize=nlmn_type) 2361 do iatom=1,natom 2362 pawrhoij(iatom)%nrhoijsel=nsel44(1,iatom) 2363 end do 2364 bsize=sum(nsel44) 2365 LIBPAW_ALLOCATE(ibuffer,(bsize)) 2366 LIBPAW_ALLOCATE(buffer,(bsize)) 2367 if (iomode == fort_binary) then 2368 read(unitfi ) ibuffer(:),buffer(:) 2369 else if (iomode == fort_formatted) then 2370 read(unitfi,*) ibuffer(:),buffer(:) 2371 end if 2372 ii=0 2373 do iatom=1,natom 2374 nselect=nsel44(1,iatom) 2375 pawrhoij(iatom)%rhoijselect(:)=0 2376 pawrhoij(iatom)%rhoijselect(1:nselect)=ibuffer(ii+1:ii+nselect) 2377 do ispden=1,nspden_in 2378 pawrhoij(iatom)%rhoijp(1:nselect,ispden)=buffer(ii+1:ii+nselect) 2379 ii=ii+nselect 2380 end do 2381 end do 2382 LIBPAW_DEALLOCATE(ibuffer) 2383 LIBPAW_DEALLOCATE(buffer) 2384 LIBPAW_DEALLOCATE(nsel44) 2385 else if (headform>=56) then 2386 LIBPAW_ALLOCATE(nsel56,(natom)) 2387 my_cplex=1;my_nspden=1;my_qphase=1 2388 if (headform==56) then 2389 if (iomode == fort_binary) then 2390 read(unitfi ) (nsel56(iatom),iatom=1,natom),my_cplex 2391 else if (iomode == fort_formatted) then 2392 read(unitfi,*) (nsel56(iatom),iatom=1,natom),my_cplex 2393 #ifdef LIBPAW_HAVE_NETCDF 2394 else if (iomode == netcdf_io) then 2395 NCF_CHECK(nf90_inq_dimid(ncid, "pawrhoij_cplex", cplex_id)) 2396 NCF_CHECK(nf90_inquire_dimension(ncid, cplex_id, len=my_cplex)) 2397 2398 NCF_CHECK(nf90_inq_varid(ncid, "rhoijsel_atoms", nsel56_id)) 2399 NCF_CHECK(nf90_get_var(ncid, nsel56_id, nsel56)) 2400 #endif 2401 end if 2402 else 2403 if (iomode == fort_binary) then 2404 read(unitfi,err=10,end=10) (nsel56(iatom),iatom=1,natom),my_cplex,my_nspden,my_qphase 2405 10 continue 2406 else if (iomode == fort_formatted) then 2407 read(unitfi,fmt=*,err=11,end=11) (nsel56(iatom),iatom=1,natom),my_cplex,my_nspden,my_qphase 2408 11 continue 2409 #ifdef LIBPAW_HAVE_NETCDF 2410 else if (iomode == netcdf_io) then 2411 NCF_CHECK(nf90_inq_dimid(ncid, "pawrhoij_cplex", cplex_id)) 2412 NCF_CHECK(nf90_inquire_dimension(ncid, cplex_id, len=my_cplex)) 2413 NCF_CHECK(nf90_inq_dimid(ncid, "pawrhoij_nspden", nspden_id)) 2414 NCF_CHECK(nf90_inquire_dimension(ncid, nspden_id, len=my_nspden)) 2415 NCF_CHECK(nf90_inq_varid(ncid, "nrhoijsel_atoms", nsel56_id)) 2416 NCF_CHECK(nf90_get_var(ncid, nsel56_id, nsel56)) 2417 if (nf90_inq_dimid(ncid, "pawrhoij_qphase", qphase_id)==NF90_NOERR) then 2418 NCF_CHECK(nf90_inquire_dimension(ncid, qphase_id, len=my_qphase)) 2419 else 2420 my_qphase=1 2421 end if 2422 #endif 2423 end if 2424 end if 2425 call pawrhoij_alloc(pawrhoij,my_cplex,my_nspden,nspinor_in,nsppol_in,typat,& 2426 & lmnsize=nlmn_type,qphase=my_qphase) 2427 do iatom=1,natom 2428 pawrhoij(iatom)%nrhoijsel=nsel56(iatom) 2429 end do 2430 bsize=sum(nsel56) 2431 LIBPAW_ALLOCATE(ibuffer,(bsize)) 2432 LIBPAW_ALLOCATE(buffer,(bsize*my_nspden*my_cplex*my_qphase)) 2433 if (iomode == fort_binary) then 2434 read(unitfi ) ibuffer(:),buffer(:) 2435 else if (iomode == fort_formatted) then 2436 read(unitfi,*) ibuffer(:),buffer(:) 2437 #ifdef LIBPAW_HAVE_NETCDF 2438 else if (iomode == netcdf_io) then 2439 if (bsize > 0) then 2440 NCF_CHECK(nf90_inq_varid(ncid, "rhoijselect_atoms", ibuffer_id)) 2441 NCF_CHECK(nf90_get_var(ncid, ibuffer_id, ibuffer)) 2442 NCF_CHECK(nf90_inq_varid(ncid, "rhoijp_atoms", buffer_id)) 2443 NCF_CHECK(nf90_get_var(ncid, buffer_id, buffer)) 2444 end if 2445 #endif 2446 end if 2447 ii=0;jj=0 2448 do iatom=1,natom 2449 nselect=nsel56(iatom) 2450 pawrhoij(iatom)%rhoijselect(:)=0 2451 pawrhoij(iatom)%rhoijselect(1:nselect)=ibuffer(ii+1:ii+nselect) 2452 ii=ii+nselect 2453 do ispden=1,my_nspden 2454 pawrhoij(iatom)%rhoijp(1:my_cplex*my_qphase*nselect,ispden)= & 2455 & buffer(jj+1:jj+my_cplex*my_qphase*nselect) 2456 jj=jj+my_cplex*my_qphase*nselect 2457 end do 2458 end do 2459 LIBPAW_DEALLOCATE(ibuffer) 2460 LIBPAW_DEALLOCATE(buffer) 2461 LIBPAW_DEALLOCATE(nsel56) 2462 end if 2463 2464 case ("W","w") ! Writing the Rhoij tab (latest format is used) 2465 2466 LIBPAW_ALLOCATE(nsel56,(natom)) 2467 my_cplex =pawrhoij(1)%cplex_rhoij 2468 my_nspden=pawrhoij(1)%nspden 2469 my_qphase=pawrhoij(1)%qphase 2470 do iatom=1,natom 2471 nsel56(iatom)=pawrhoij(iatom)%nrhoijsel 2472 end do 2473 bsize=sum(nsel56) 2474 2475 if (iomode == fort_binary) then 2476 write(unitfi ) (nsel56(iatom),iatom=1,natom),my_cplex,my_nspden,my_qphase 2477 else if (iomode == fort_formatted) then 2478 write(unitfi,*) (nsel56(iatom),iatom=1,natom),my_cplex,my_nspden,my_qphase 2479 #ifdef LIBPAW_HAVE_NETCDF 2480 else if (iomode == netcdf_io) then 2481 ncerr = nf90_redef(ncid) 2482 2483 ! Define dimensions. 2484 ncerr = nf90_inq_dimid(ncid, "number_of_atoms", natom_id) 2485 if (ncerr /= nf90_noerr) then 2486 NCF_CHECK(nf90_def_dim(ncid, "number_of_atoms", natom, natom_id)) 2487 end if 2488 ncerr = nf90_inq_varid(ncid, "nrhoijsel_atoms", nsel56_id) 2489 if (ncerr /= nf90_noerr) then 2490 NCF_CHECK(nf90_def_var(ncid, "nrhoijsel_atoms", NF90_INT, natom_id, nsel56_id)) 2491 end if 2492 ncerr = nf90_inq_dimid(ncid, "pawrhoij_cplex", cplex_id) 2493 if (ncerr /= nf90_noerr) then 2494 NCF_CHECK(nf90_def_dim(ncid, "pawrhoij_cplex", my_cplex, cplex_id)) 2495 end if 2496 ncerr = nf90_inq_dimid(ncid, "pawrhoij_nspden", nspden_id) 2497 if (ncerr /= nf90_noerr) then 2498 NCF_CHECK(nf90_def_dim(ncid, "pawrhoij_nspden", my_nspden, nspden_id)) 2499 end if 2500 ncerr = nf90_inq_dimid(ncid, "pawrhoij_qphase", qphase_id) 2501 if (ncerr /= nf90_noerr) then 2502 NCF_CHECK(nf90_def_dim(ncid, "pawrhoij_qphase", my_qphase, qphase_id)) 2503 end if 2504 if (bsize > 0) then 2505 ncerr = nf90_inq_dimid(ncid, "rhoijselect_atoms_dim", bsize_id) 2506 if (ncerr /= nf90_noerr) then 2507 NCF_CHECK(nf90_def_dim(ncid, "rhoijselect_atoms_dim", bsize, bsize_id)) 2508 end if 2509 ncerr = nf90_inq_dimid(ncid, "rhoijp_atoms_dim", bufsize_id) 2510 if (ncerr /= nf90_noerr) then 2511 NCF_CHECK(nf90_def_dim(ncid, "rhoijp_atoms_dim", bsize*my_nspden*my_qphase*my_cplex, bufsize_id)) 2512 end if 2513 ncerr = nf90_inq_varid(ncid, "rhoijselect_atoms", ibuffer_id) 2514 if (ncerr /= nf90_noerr) then 2515 NCF_CHECK(nf90_def_var(ncid, "rhoijselect_atoms", NF90_INT, bsize_id, ibuffer_id)) 2516 end if 2517 ncerr = nf90_inq_varid(ncid, "rhoijp_atoms", buffer_id) 2518 if (ncerr /= nf90_noerr) then 2519 NCF_CHECK(nf90_def_var(ncid, "rhoijp_atoms", NF90_DOUBLE, bufsize_id, buffer_id)) 2520 end if 2521 else 2522 ! This happens in v5[40] and bsize == 0 corresponds to NC_UNLIMITED 2523 LIBPAW_COMMENT("All rhoij entries are zero. No netcdf entry produced") 2524 end if 2525 2526 ! Write nsel56 2527 NCF_CHECK(nf90_enddef(ncid)) 2528 NCF_CHECK(nf90_put_var(ncid, nsel56_id, nsel56)) 2529 #endif 2530 end if 2531 2532 LIBPAW_ALLOCATE(ibuffer,(bsize)) 2533 LIBPAW_ALLOCATE(buffer,(bsize*my_nspden*my_cplex*my_qphase)) 2534 ii=0;jj=0 2535 do iatom=1,natom 2536 nselect=nsel56(iatom) 2537 ibuffer(ii+1:ii+nselect)=pawrhoij(iatom)%rhoijselect(1:nselect) 2538 ii=ii+nselect 2539 do ispden=1,my_nspden 2540 buffer(jj+1:jj+my_cplex*my_qphase*nselect)= & 2541 & pawrhoij(iatom)%rhoijp(1:my_cplex*my_qphase*nselect,ispden) 2542 jj=jj+my_cplex*my_qphase*nselect 2543 end do 2544 end do 2545 if (iomode == fort_binary) then 2546 write(unitfi ) ibuffer(:),buffer(:) 2547 else if (iomode == fort_formatted) then 2548 write(unitfi,*) ibuffer(:),buffer(:) 2549 #ifdef LIBPAW_HAVE_NETCDF 2550 else if (iomode == netcdf_io) then 2551 if (bsize > 0) then 2552 NCF_CHECK(nf90_put_var(ncid, ibuffer_id, ibuffer)) 2553 NCF_CHECK(nf90_put_var(ncid, buffer_id, buffer)) 2554 end if 2555 #endif 2556 end if 2557 LIBPAW_DEALLOCATE(ibuffer) 2558 LIBPAW_DEALLOCATE(buffer) 2559 LIBPAW_DEALLOCATE(nsel56) 2560 2561 case ("E","e") ! Echoing the Rhoij tab 2562 2563 my_natinc=1; if(natom>1) my_natinc=natom-1 2564 my_qphase=pawrhoij(1)%qphase 2565 nselect=maxval(pawrhoij(:)%nrhoijsel) 2566 if (PRESENT(natinc)) my_natinc = natinc ! user-defined increment. 2567 LIBPAW_ALLOCATE(ibuffer,(0)) 2568 nselect=maxval(pawrhoij(:)%nrhoijsel) 2569 if (my_qphase==2) then 2570 LIBPAW_POINTER_ALLOCATE(rhoij_tmp,(2*nselect)) 2571 end if 2572 do iatom=1,my_natom,my_natinc 2573 iatom_tot=iatom;if(paral_atom)iatom_tot=my_atmtab(iatom) 2574 my_cplex =pawrhoij(iatom)%cplex_rhoij 2575 my_nspden=pawrhoij(iatom)%nspden 2576 my_qphase=pawrhoij(iatom)%qphase 2577 nselect=pawrhoij(iatom)%nrhoijsel 2578 do ispden=1,pawrhoij(iatom)%nspden 2579 if (my_qphase==1) then 2580 my_cplex_eff=my_cplex 2581 rhoij_tmp => pawrhoij(iatom)%rhoijp(1:my_cplex*nselect,ispden) 2582 else 2583 my_cplex_eff=2 2584 rhoij_tmp=zero 2585 jj=my_cplex*pawrhoij(iatom)%lmn2_size 2586 if (my_cplex==1) then 2587 do ii=1,nselect 2588 rhoij_tmp(2*ii-1)=pawrhoij(iatom)%rhoijp(ii,ispden) 2589 rhoij_tmp(2*ii )=pawrhoij(iatom)%rhoijp(jj+ii,ispden) 2590 end do 2591 else 2592 do ii=1,nselect 2593 rhoij_tmp(2*ii-1)=pawrhoij(iatom)%rhoijp(2*ii-1,ispden) & 2594 & -pawrhoij(iatom)%rhoijp(jj+2*ii ,ispden) 2595 rhoij_tmp(2*ii )=pawrhoij(iatom)%rhoijp(2*ii ,ispden) & 2596 & +pawrhoij(iatom)%rhoijp(jj+2*ii-1,ispden) 2597 end do 2598 end if 2599 end if 2600 write(unitfi, '(a,i4,a,i1,a)' ) ' rhoij(',iatom_tot,',',ispden,')= (max 12 non-zero components will be written)' 2601 call pawio_print_ij(unitfi,rhoij_tmp,nselect,my_cplex_eff,& 2602 & pawrhoij(iatom)%lmn_size,-1,ibuffer,1,0,& 2603 & pawrhoij(iatom)%rhoijselect,-1.d0,1,& 2604 & opt_sym=2,mode_paral='PERS') 2605 end do ! end nspden do 2606 end do ! end iatom do 2607 LIBPAW_DEALLOCATE(ibuffer) 2608 if (my_qphase==2) then 2609 LIBPAW_POINTER_DEALLOCATE(rhoij_tmp) 2610 end if 2611 2612 case ("D","d") ! Debug 2613 2614 write(unitfi,'(a,i4)' ) 'size pawmkrhoij , natom = ' , my_natom 2615 my_natinc=1; if(natom>1) my_natinc=natom-1 2616 if (PRESENT(natinc)) my_natinc = natinc ! user-defined increment. 2617 LIBPAW_ALLOCATE(ibuffer,(0)) 2618 do iatom=1,my_natom,my_natinc 2619 iatom_tot=iatom;if(paral_atom) iatom_tot=my_atmtab(iatom) 2620 if (iatom_tot/=1) cycle 2621 write(unitfi,'(a,i4,a)' ) ' ******* rhoij (Atom # ',iatom_tot,' ********)' 2622 write(unitfi,'(a,i4,a,i4)' ) 'cplex_rhoij=',pawrhoij(iatom)%cplex_rhoij, ' nselect=', pawrhoij(iatom)%nrhoijsel 2623 write(unitfi,'(a,i4,a,i4)' ) 'nspden=',pawrhoij(iatom)%nspden, ' lmn2size=',pawrhoij(iatom)%lmn2_size 2624 write(unitfi,'(a,i4,a,i4)' ) 'lmnmix=',pawrhoij(iatom)%lmnmix_sz, ' ngrhoij=',pawrhoij(iatom)%ngrhoij 2625 write(unitfi,'(a,i4,a,i4)' ) 'use_rhoijres=',pawrhoij(iatom)%use_rhoijres, & 2626 & 'use_rhoij_=',pawrhoij(iatom)%use_rhoij_ 2627 write(unitfi,'(a,i4,a,i4)' ) 'itypat=',pawrhoij(iatom)%itypat, ' lmn_size=',pawrhoij(iatom)%lmn_size 2628 write(unitfi,'(a,i4,a,i4)' ) 'nsppol=',pawrhoij(iatom)%nsppol, ' nspinor=',pawrhoij(iatom)%nspinor 2629 write(unitfi,'(a,i4)' ) 'qphase=',pawrhoij(iatom)%qphase 2630 cplex=pawrhoij(iatom)%cplex_rhoij 2631 qphase=pawrhoij(iatom)%qphase 2632 lmn2_size=pawrhoij(iatom)%lmn2_size 2633 do i2=1,pawrhoij(iatom)%nrhoijsel 2634 write(unitfi,'(a,i4,a,i4,a,i4,a,f9.5)') 'rhoijselect(,',i2,')=',& 2635 & pawrhoij(iatom)%rhoijselect(i2) 2636 end do 2637 if (pawrhoij(iatom)%ngrhoij>0) then 2638 ngrhoijmx=2; if (pawrhoij(iatom)%ngrhoij<ngrhoijmx) ngrhoijmx=pawrhoij(iatom)%ngrhoij 2639 do ispden=1,pawrhoij(iatom)%nspden 2640 do i1=ngrhoijmx,ngrhoijmx 2641 do ii=1,qphase 2642 do i2=(ii-1)*cplex*lmn2_size+cplex*lmn2_size,(ii-1)*cplex*lmn2_size+cplex*lmn2_size 2643 write(unitfi,'(a,i4,a,i4,a,i4,a,f9.5)') ' grhoij(,',i1,',',i2,',',ispden,')=',& 2644 & pawrhoij(iatom)%grhoij(i1,i2,ispden) 2645 end do 2646 end do 2647 end do 2648 end do 2649 call libpaw_flush(unitfi) 2650 end if 2651 if (pawrhoij(iatom)%use_rhoijres>0) then 2652 do ispden=1,pawrhoij(iatom)%nspden 2653 do ii=1,qphase 2654 do i2=(ii-1)*cplex*lmn2_size+cplex*lmn2_size,(ii-1)*cplex*lmn2_size+cplex*lmn2_size 2655 write(unitfi,'(a,i4,a,i4,a,f9.5)') ' rhoijres(,',i2,',ispden=',ispden,')=',& 2656 & pawrhoij(iatom)%rhoijres(i2,ispden) 2657 end do 2658 end do 2659 end do 2660 call libpaw_flush(unitfi) 2661 end if 2662 if (pawrhoij(iatom)%nrhoijsel>0) then 2663 do ispden=1,pawrhoij(iatom)%nspden 2664 do ii=1,qphase 2665 do i2=(ii-1)*cplex*lmn2_size+cplex*pawrhoij(iatom)%nrhoijsel, & 2666 & (ii-1)*cplex*lmn2_size+cplex*pawrhoij(iatom)%nrhoijsel 2667 write(unitfi,'(a,i4,a,i4,a,f9.5)') ' rhoijp!(nrhoijselec=,',i2,',ispden=',ispden,')=',& 2668 & pawrhoij(iatom)%rhoijp(i2,ispden) 2669 end do 2670 end do 2671 end do 2672 call libpaw_flush(unitfi) 2673 end if 2674 if (pawrhoij(iatom)%use_rhoij_>0) then 2675 size_rhoij2=size(pawrhoij(iatom)%rhoij_,2) 2676 do ispden=1,size_rhoij2 2677 do ii=1,qphase 2678 do i2=(ii-1)*cplex*lmn2_size+cplex*lmn2_size,(ii-1)*cplex*lmn2_size+cplex*lmn2_size 2679 write(unitfi,'(a,i4,a,i4,a,f9.5)') ' rhoij_(,',i2,',ispden=',ispden,')=',& 2680 & pawrhoij(iatom)%rhoij_(i2,ispden) 2681 end do 2682 end do 2683 end do 2684 end if 2685 call libpaw_flush(unitfi) 2686 if (pawrhoij(iatom)%lmnmix_sz>0) then 2687 write(unitfi,'(a)') 'kpawmix ' 2688 write(unitfi,'(i4,i4,i4,i4)') pawrhoij(iatom)%kpawmix(:) 2689 end if 2690 call libpaw_flush(unitfi) 2691 end do 2692 2693 case default 2694 msg='Wrong rdwr_mode'//TRIM(rdwr_mode) 2695 LIBPAW_ERROR(msg) 2696 2697 end select 2698 2699 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
4471 subroutine pawrhoij_isendreceive_fillbuffer(pawrhoij,atmtab_send, atm_indx_send,nrhoij_send,& 4472 & buf_int,buf_int_size,buf_dp,buf_dp_size) 4473 4474 !Arguments ------------------------------------ 4475 !scalars 4476 integer,intent(out) :: buf_dp_size,buf_int_size 4477 integer,intent(in) :: nrhoij_send 4478 !arrays 4479 integer,intent(in) :: atmtab_send(:),atm_indx_send(:) 4480 integer,intent(out),allocatable :: buf_int(:) 4481 real(dp),intent(out),allocatable :: buf_dp(:) 4482 type(pawrhoij_type),target,intent(in) :: pawrhoij(:) 4483 4484 !Local variables------------------------------- 4485 !scalars 4486 integer :: cplex,ii,indx_int,indx_dp, iatom_tot,irhoij,irhoij_send,isp,jj,lmn2_size,lmnmix 4487 integer :: ngrhoij,nselect,nspden,qphase,rhoij_size2 4488 integer :: use_rhoijp,use_rhoijres,use_rhoij_ 4489 character(len=500) :: msg 4490 type(pawrhoij_type),pointer :: pawrhoij1 4491 !arrays 4492 4493 ! ********************************************************************* 4494 4495 !Compute sizes of buffers 4496 buf_int_size=0;buf_dp_size=0 4497 nselect=0;lmnmix=0;ngrhoij=0;rhoij_size2=0 4498 use_rhoijp=0;use_rhoijres=0;use_rhoij_=0 4499 do irhoij_send=1,nrhoij_send 4500 iatom_tot=atmtab_send(irhoij_send) 4501 irhoij=atm_indx_send(iatom_tot) 4502 if (irhoij == -1) then 4503 msg="Error in pawrhoij_isendreceive_fillbuffer atom not found" 4504 LIBPAW_BUG(msg) 4505 end if 4506 pawrhoij1=>pawrhoij(irhoij) 4507 cplex =pawrhoij1%cplex_rhoij 4508 qphase =pawrhoij1%qphase 4509 lmn2_size=pawrhoij1%lmn2_size 4510 nspden =pawrhoij1%nspden 4511 lmnmix=pawrhoij1%lmnmix_sz 4512 ngrhoij=pawrhoij1%ngrhoij 4513 use_rhoijp=pawrhoij1%use_rhoijp 4514 use_rhoijres=pawrhoij1%use_rhoijres 4515 use_rhoij_=pawrhoij1%use_rhoij_ 4516 buf_int_size=buf_int_size+16 4517 if (use_rhoijp>0) then 4518 nselect=pawrhoij1%nrhoijsel 4519 buf_int_size=buf_int_size+nselect 4520 buf_dp_size=buf_dp_size + cplex*qphase*nselect*nspden 4521 end if 4522 if (lmnmix>0) buf_int_size=buf_int_size+lmnmix 4523 if (ngrhoij>0) buf_dp_size=buf_dp_size + cplex*qphase*lmn2_size*nspden*ngrhoij 4524 if (use_rhoijres>0) buf_dp_size=buf_dp_size + cplex*qphase*lmn2_size*nspden 4525 if (use_rhoij_>0) then 4526 rhoij_size2=size(pawrhoij1%rhoij_,dim=2) 4527 buf_dp_size=buf_dp_size + cplex*qphase*lmn2_size*rhoij_size2 4528 end if 4529 end do 4530 4531 !Fill input buffers 4532 LIBPAW_ALLOCATE(buf_int,(buf_int_size)) 4533 LIBPAW_ALLOCATE(buf_dp,(buf_dp_size)) 4534 indx_int=1;indx_dp =1 4535 lmnmix=0;ngrhoij=0;nselect=0;rhoij_size2=0 4536 use_rhoijp=0;use_rhoijres=0;use_rhoij_=0 4537 do irhoij_send=1,nrhoij_send 4538 iatom_tot=atmtab_send(irhoij_send) 4539 irhoij=atm_indx_send(iatom_tot) 4540 pawrhoij1=>pawrhoij(irhoij) 4541 cplex =pawrhoij1%cplex_rhoij 4542 qphase =pawrhoij1%qphase 4543 lmn2_size=pawrhoij1%lmn2_size 4544 nspden =pawrhoij1%nspden 4545 lmnmix=pawrhoij1%lmnmix_sz 4546 ngrhoij=pawrhoij1%ngrhoij 4547 use_rhoijp=pawrhoij1%use_rhoijp 4548 nselect=pawrhoij1%nrhoijsel 4549 use_rhoijres=pawrhoij1%use_rhoijres 4550 use_rhoij_ =pawrhoij1%use_rhoij_ 4551 if (use_rhoij_ > 0) then 4552 rhoij_size2 =size(pawrhoij1%rhoij_,dim=2) 4553 end if 4554 buf_int(indx_int)=atmtab_send(irhoij_send) ;indx_int=indx_int+1 4555 buf_int(indx_int)=cplex ;indx_int=indx_int+1 4556 buf_int(indx_int)=qphase ;indx_int=indx_int+1 4557 buf_int(indx_int)=lmn2_size ;indx_int=indx_int+1 4558 buf_int(indx_int)=nspden ;indx_int=indx_int+1 4559 buf_int(indx_int)=nselect ;indx_int=indx_int+1 4560 buf_int(indx_int)=lmnmix ;indx_int=indx_int+1 4561 buf_int(indx_int)=ngrhoij ;indx_int=indx_int+1 4562 buf_int(indx_int)=use_rhoijp ;indx_int=indx_int+1 4563 buf_int(indx_int)=use_rhoijres ;indx_int=indx_int+1 4564 buf_int(indx_int)=use_rhoij_ ;indx_int=indx_int+1 4565 buf_int(indx_int)=rhoij_size2 ;indx_int=indx_int+1 4566 buf_int(indx_int)=pawrhoij1%itypat ;indx_int=indx_int+1 4567 buf_int(indx_int)=pawrhoij1%lmn_size ;indx_int=indx_int+1 4568 buf_int(indx_int)=pawrhoij1%nsppol ;indx_int=indx_int+1 4569 buf_int(indx_int)=pawrhoij1%nspinor ;indx_int=indx_int+1 4570 if (use_rhoijp>0) then 4571 buf_int(indx_int:indx_int+nselect-1)=pawrhoij1%rhoijselect(1:nselect) 4572 indx_int=indx_int+nselect 4573 do isp=1,nspden 4574 do ii=1,qphase 4575 jj=(ii-1)*cplex*lmn2_size 4576 buf_dp(indx_dp:indx_dp+cplex*nselect-1)=pawrhoij1%rhoijp(jj+1:jj+cplex*nselect,isp) 4577 indx_dp=indx_dp+cplex*nselect 4578 end do 4579 end do 4580 end if 4581 if (lmnmix>0) then 4582 buf_int(indx_int:indx_int+lmnmix-1)=pawrhoij1%kpawmix(1:lmnmix) 4583 indx_int=indx_int+lmnmix 4584 end if 4585 if (ngrhoij>0) then 4586 do isp=1,nspden 4587 do ii=1,cplex*qphase*lmn2_size 4588 buf_dp(indx_dp:indx_dp+ngrhoij-1)=pawrhoij1%grhoij(1:ngrhoij,ii,isp) 4589 indx_dp=indx_dp+ngrhoij 4590 end do 4591 end do 4592 end if 4593 if (use_rhoijres>0) then 4594 do isp=1,nspden 4595 buf_dp(indx_dp:indx_dp+cplex*qphase*lmn2_size-1)=pawrhoij1%rhoijres(1:cplex*qphase*lmn2_size,isp) 4596 indx_dp=indx_dp+cplex*qphase*lmn2_size 4597 end do 4598 end if 4599 if (use_rhoij_>0) then 4600 do isp=1,rhoij_size2 4601 buf_dp(indx_dp:indx_dp+cplex*qphase*lmn2_size-1)=pawrhoij1%rhoij_(1:cplex*qphase*lmn2_size,isp) 4602 indx_dp=indx_dp+cplex*qphase*lmn2_size 4603 end do 4604 end if 4605 end do !irhoij_send 4606 4607 !Check 4608 if ((indx_int-1/=buf_int_size).or.(indx_dp-1/=buf_dp_size)) then 4609 write(msg,'(a,i10,a,i10)') 'Wrong buffer sizes: buf_int_size=',buf_int_size,' buf_dp_size=',buf_dp_size 4610 LIBPAW_BUG(msg) 4611 end if 4612 4613 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
4332 subroutine pawrhoij_isendreceive_getbuffer(pawrhoij,nrhoij_send,atm_indx_recv,buf_int,buf_dp) 4333 4334 !Arguments ------------------------------------ 4335 !scalars 4336 integer,intent(in) :: nrhoij_send 4337 !arrays 4338 integer,intent(in) ::atm_indx_recv(:),buf_int(:) 4339 real(dp),intent(in) :: buf_dp(:) 4340 type(pawrhoij_type),target,intent(inout) :: pawrhoij(:) 4341 4342 !Local variables------------------------------- 4343 !scalars 4344 integer :: buf_dp_size,buf_int_size,cplex,ii,indx_int,indx_dp,iatom_tot,irhoij_send 4345 integer :: isp,jj,jrhoij,lmn2_size,lmnmix,ngrhoij,nselect,nspden,qphase,rhoij_size2 4346 integer :: use_rhoijp,use_rhoijres,use_rhoij_ 4347 character(len=500) :: msg 4348 type(pawrhoij_type),pointer :: pawrhoij1 4349 !arrays 4350 4351 ! ********************************************************************* 4352 4353 buf_int_size=size(buf_int) 4354 buf_dp_size=size(buf_dp) 4355 indx_int=1;indx_dp=1 4356 4357 do irhoij_send=1,nrhoij_send 4358 4359 iatom_tot=buf_int(indx_int) ; indx_int=indx_int+1 4360 jrhoij= atm_indx_recv(iatom_tot) 4361 if (jrhoij==-1) then 4362 msg="Error in pawrhoij_isendreceive_getbuffer atom not found" 4363 LIBPAW_BUG(msg) 4364 end if 4365 pawrhoij1=>pawrhoij(jrhoij) 4366 4367 cplex =buf_int(indx_int) ;indx_int=indx_int+1 4368 qphase =buf_int(indx_int) ;indx_int=indx_int+1 4369 lmn2_size =buf_int(indx_int) ;indx_int=indx_int+1 4370 nspden =buf_int(indx_int) ;indx_int=indx_int+1 4371 nselect =buf_int(indx_int) ;indx_int=indx_int+1 4372 lmnmix =buf_int(indx_int) ;indx_int=indx_int+1 4373 ngrhoij =buf_int(indx_int) ;indx_int=indx_int+1 4374 use_rhoijp =buf_int(indx_int) ;indx_int=indx_int+1 4375 use_rhoijres=buf_int(indx_int) ;indx_int=indx_int+1 4376 use_rhoij_ =buf_int(indx_int) ;indx_int=indx_int+1 4377 rhoij_size2 =buf_int(indx_int) ;indx_int=indx_int+1 4378 pawrhoij1%itypat=buf_int(indx_int) ;indx_int=indx_int+1 4379 pawrhoij1%lmn_size=buf_int(indx_int) ;indx_int=indx_int+1 4380 pawrhoij1%nsppol=buf_int(indx_int) ;indx_int=indx_int+1 4381 pawrhoij1%nspinor=buf_int(indx_int) ;indx_int=indx_int+1 4382 pawrhoij1%cplex_rhoij=cplex 4383 pawrhoij1%qphase=qphase 4384 pawrhoij1%lmn2_size=lmn2_size 4385 pawrhoij1%nspden=nspden 4386 pawrhoij1%nrhoijsel=nselect 4387 pawrhoij1%lmnmix_sz=lmnmix 4388 pawrhoij1%ngrhoij=ngrhoij 4389 pawrhoij1%use_rhoijp=use_rhoijp 4390 pawrhoij1%use_rhoijres=use_rhoijres 4391 pawrhoij1%use_rhoij_=use_rhoij_ 4392 if (use_rhoijp>0) then 4393 LIBPAW_ALLOCATE(pawrhoij1%rhoijselect,(lmn2_size)) 4394 pawrhoij1%rhoijselect(1:nselect)=buf_int(indx_int:indx_int+nselect-1) 4395 if (nselect < lmn2_size )pawrhoij1%rhoijselect(nselect+1:lmn2_size)=zero 4396 indx_int=indx_int+nselect 4397 LIBPAW_ALLOCATE(pawrhoij1%rhoijp,(cplex*qphase*lmn2_size,nspden)) 4398 do isp=1,nspden 4399 do ii=1,qphase 4400 jj=(ii-1)*cplex*lmn2_size 4401 pawrhoij1%rhoijp(jj+1:jj+cplex*nselect,isp)=buf_dp(indx_dp:indx_dp+cplex*nselect-1) 4402 if (nselect<lmn2_size)pawrhoij1%rhoijp(jj+cplex*nselect+1:jj+cplex*lmn2_size,isp)=zero 4403 indx_dp=indx_dp+cplex*nselect 4404 end do 4405 end do 4406 end if 4407 if (lmnmix>0) then 4408 LIBPAW_ALLOCATE(pawrhoij1%kpawmix,(lmnmix)) 4409 pawrhoij1%kpawmix(1:lmnmix)=buf_int(indx_int:indx_int+lmnmix-1) 4410 indx_int=indx_int+lmnmix 4411 end if 4412 if (ngrhoij>0) then 4413 LIBPAW_ALLOCATE(pawrhoij1%grhoij,(ngrhoij,cplex*qphase*lmn2_size,nspden)) 4414 do isp=1,nspden 4415 do ii=1,cplex*qphase*lmn2_size 4416 pawrhoij1%grhoij(1:ngrhoij,ii,isp)=buf_dp(indx_dp:indx_dp+ngrhoij-1) 4417 indx_dp=indx_dp+ngrhoij 4418 end do 4419 end do 4420 end if 4421 if (use_rhoijres>0) then 4422 LIBPAW_ALLOCATE(pawrhoij1%rhoijres,(cplex*qphase*lmn2_size,nspden)) 4423 do isp=1,nspden 4424 pawrhoij1%rhoijres(1:cplex*qphase*lmn2_size,isp)=buf_dp(indx_dp:indx_dp+cplex*qphase*lmn2_size-1) 4425 indx_dp=indx_dp+cplex*qphase*lmn2_size 4426 end do 4427 end if 4428 if (use_rhoij_>0) then 4429 LIBPAW_ALLOCATE(pawrhoij1%rhoij_,(cplex*qphase*lmn2_size,rhoij_size2)) 4430 do isp=1,rhoij_size2 4431 pawrhoij1%rhoij_(1:cplex*qphase*lmn2_size,isp)=buf_dp(indx_dp:indx_dp+cplex*qphase*lmn2_size-1) 4432 indx_dp=indx_dp+cplex*qphase*lmn2_size 4433 end do 4434 end if 4435 end do !irhoij_send 4436 if ((indx_int/=1+buf_int_size).or.(indx_dp/=1+buf_dp_size)) then 4437 write(msg,'(a,i10,a,i10)') 'Wrong buffer sizes: buf_int_size=',buf_int_size,' buf_dp_size=',buf_dp_size 4438 LIBPAW_BUG(msg) 4439 end if 4440 4441 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
2881 subroutine pawrhoij_mpisum_unpacked_1D(pawrhoij,comm1,comm2) 2882 2883 !Arguments --------------------------------------------- 2884 !scalars 2885 integer,intent(in) :: comm1 2886 integer,optional,intent(in) :: comm2 2887 !arrays 2888 type(pawrhoij_type),intent(inout) :: pawrhoij(:) 2889 2890 !Local variables --------------------------------------- 2891 !scalars 2892 integer :: bufdim,iatom,ierr,isppol,jdim,nsp2,natom 2893 integer :: nproc1,nproc2 2894 !character(len=500) :: msg 2895 !arrays 2896 integer,allocatable :: dimlmn(:) 2897 real(dp),allocatable :: buffer(:) 2898 2899 !************************************************************************ 2900 2901 natom=SIZE(pawrhoij);if (natom==0) return 2902 2903 nproc1 = xmpi_comm_size(comm1) 2904 nproc2=1; if (PRESENT(comm2)) nproc2 = xmpi_comm_size(comm2) 2905 if (nproc1==1.and.nproc2==1) RETURN 2906 2907 !Fill the MPI buffer from the local rhoij_ 2908 LIBPAW_ALLOCATE(dimlmn,(natom)) 2909 dimlmn(1:natom)=pawrhoij(1:natom)%cplex_rhoij & 2910 & *pawrhoij(1:natom)%qphase & 2911 & *pawrhoij(1:natom)%lmn2_size 2912 nsp2=pawrhoij(1)%nsppol; if (pawrhoij(1)%nspden==4) nsp2=4 2913 bufdim=sum(dimlmn)*nsp2 2914 LIBPAW_ALLOCATE(buffer,(bufdim)) 2915 jdim=0 2916 do iatom=1,natom 2917 do isppol=1,nsp2 2918 buffer(jdim+1:jdim+dimlmn(iatom))=pawrhoij(iatom)%rhoij_(:,isppol) 2919 jdim=jdim+dimlmn(iatom) 2920 end do 2921 end do 2922 2923 !Build sum of pawrhoij%rhoij_ 2924 call xmpi_sum(buffer,comm1,ierr) ! Sum over the first communicator. 2925 if (PRESENT(comm2)) then 2926 call xmpi_sum(buffer,comm2,ierr) ! Sum over the second communicator. 2927 end if 2928 2929 !Unpack the MPI packet filling rhoij_ 2930 jdim=0 2931 do iatom=1,natom 2932 do isppol=1,nsp2 2933 pawrhoij(iatom)%rhoij_(:,isppol)=buffer(jdim+1:jdim+dimlmn(iatom)) 2934 jdim=jdim+dimlmn(iatom) 2935 end do 2936 end do 2937 2938 LIBPAW_DEALLOCATE(buffer) 2939 LIBPAW_DEALLOCATE(dimlmn) 2940 2941 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
2967 subroutine pawrhoij_mpisum_unpacked_2D(pawrhoij,comm1,comm2) 2968 2969 !Arguments --------------------------------------------- 2970 !scalars 2971 integer,intent(in) :: comm1 2972 integer,optional,intent(in) :: comm2 2973 !arrays 2974 type(pawrhoij_type),intent(inout) :: pawrhoij(:,:) 2975 2976 !Local variables --------------------------------------- 2977 !scalars 2978 integer :: bufdim,iatom,ierr,irhoij,isppol,jdim,nsp2,natom,nrhoij 2979 integer :: nproc1,nproc2 2980 !character(len=500) :: msg 2981 !arrays 2982 integer,allocatable :: dimlmn(:,:) 2983 real(dp),allocatable :: buffer(:) 2984 2985 !************************************************************************ 2986 2987 natom =SIZE(pawrhoij,1);if (natom ==0) return 2988 nrhoij=SIZE(pawrhoij,2);if (nrhoij==0) return 2989 2990 nproc1 = xmpi_comm_size(comm1) 2991 nproc2=1; if (PRESENT(comm2)) nproc2 = xmpi_comm_size(comm2) 2992 if (nproc1==1.and.nproc2==1) RETURN 2993 2994 !Fill the MPI buffer from the local rhoij_ 2995 LIBPAW_ALLOCATE(dimlmn,(natom,nrhoij)) 2996 dimlmn(1:natom,1:nrhoij)=pawrhoij(1:natom,1:nrhoij)%cplex_rhoij & 2997 & *pawrhoij(1:natom,1:nrhoij)%qphase & 2998 & *pawrhoij(1:natom,1:nrhoij)%lmn2_size 2999 nsp2=pawrhoij(1,1)%nsppol; if (pawrhoij(1,1)%nspden==4) nsp2=4 3000 bufdim=sum(dimlmn)*nsp2 3001 LIBPAW_ALLOCATE(buffer,(bufdim)) 3002 jdim=0 3003 do irhoij=1,nrhoij 3004 do iatom=1,natom 3005 do isppol=1,nsp2 3006 buffer(jdim+1:jdim+dimlmn(iatom,irhoij))=pawrhoij(iatom,irhoij)%rhoij_(:,isppol) 3007 jdim=jdim+dimlmn(iatom,irhoij) 3008 end do 3009 end do 3010 end do 3011 3012 !Build sum of pawrhoij%rhoij_ 3013 call xmpi_sum(buffer,comm1,ierr) ! Sum over the first communicator. 3014 if (PRESENT(comm2)) then 3015 call xmpi_sum(buffer,comm2,ierr) ! Sum over the second communicator. 3016 end if 3017 3018 !Unpack the MPI packet filling rhoij_ 3019 jdim=0 3020 do irhoij=1,nrhoij 3021 do iatom=1,natom 3022 do isppol=1,nsp2 3023 pawrhoij(iatom,irhoij)%rhoij_(:,isppol)=buffer(jdim+1:jdim+dimlmn(iatom,irhoij)) 3024 jdim=jdim+dimlmn(iatom,irhoij) 3025 end do 3026 end do 3027 end do 3028 3029 LIBPAW_DEALLOCATE(buffer) 3030 LIBPAW_DEALLOCATE(dimlmn) 3031 3032 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
3475 subroutine pawrhoij_symrhoij(pawrhoij,pawrhoij_unsym,choice,gprimd,indsym,ipert,natom,nsym,& 3476 & ntypat,optrhoij,pawang,pawprtvol,pawtab,rprimd,symafm,symrec,typat, & 3477 & mpi_atmtab,comm_atom,qphon,use_zeromag) ! optional arguments (parallelism) 3478 3479 !Arguments --------------------------------------------- 3480 !scalars 3481 integer,intent(in) :: choice,ipert,natom,nsym,ntypat,optrhoij,pawprtvol 3482 integer,optional,intent(in) :: comm_atom 3483 logical,optional,intent(in) :: use_zeromag 3484 type(pawang_type),intent(in) :: pawang 3485 !arrays 3486 integer,intent(in) :: indsym(4,nsym,natom) 3487 integer,optional,target,intent(in) :: mpi_atmtab(:) 3488 integer,intent(in) :: symafm(nsym),symrec(3,3,nsym),typat(natom) 3489 real(dp),intent(in) :: gprimd(3,3),rprimd(3,3) 3490 real(dp),intent(in),optional :: qphon(3) 3491 type(pawrhoij_type),intent(inout) :: pawrhoij(:) 3492 type(pawrhoij_type),target,intent(inout) :: pawrhoij_unsym(:) 3493 type(pawtab_type),target,intent(in) :: pawtab(ntypat) 3494 3495 !Local variables --------------------------------------- 3496 !scalars 3497 integer :: at_indx,cplex_eff,cplex_rhoij,iafm,iatm,iatom,il,il0,ilmn,iln,iln0,ilpm,indexi 3498 integer :: indexii,indexj,indexjj,indexjj0,indexk,indexkc,indexkc_q,iplex,iq,iq0 3499 integer :: irhoij,irot,ishift2,ishift3,ishift4,ispden,itypat 3500 integer :: j0lmn,jl,jl0,jlmn,jln,jln0,jlpm,jrhoij,jspden,klmn,klmn1,klmn1q,kspden 3501 integer :: lmn_size,lmn2_size,mi,mj,my_comm_atom,mu,mua,mub,mushift 3502 integer :: natinc,ngrhoij,nrhoij,nrhoij1,nrhoij_unsym 3503 integer :: nselect,nu,nushift,qphase,sz1,sz2 3504 logical,parameter :: afm_noncoll=.true. ! TRUE if antiferro symmetries are used with non-collinear magnetism 3505 logical :: use_zeromag_ 3506 real(dp) :: arg,factafm,ro,syma,zarot2 3507 logical :: antiferro,has_qphase,my_atmtab_allocated,noncoll 3508 logical :: paral_atom,paral_atom_unsym,use_afm,use_res 3509 character(len=8) :: pertstrg,wrt_mode 3510 character(len=500) :: msg 3511 !arrays 3512 integer,parameter :: alpha(6)=(/1,2,3,3,3,2/),beta(6)=(/1,2,3,2,1,1/) 3513 integer :: nsym_used(2) 3514 integer, pointer :: indlmn(:,:) 3515 integer,pointer :: my_atmtab(:) 3516 real(dp) :: fact(2),factsym(2),phase(2),rhoijc(2),rotmag(2,3,2),rotrho(2,2,2) 3517 real(dp) :: summag(2,3,2),sumrho(2,2,2),sum1(2),work1(2,3,3),xsym(3) 3518 real(dp),allocatable :: rotgr(:,:,:,:),rotmaggr(:,:,:,:),sumgr(:,:,:),summaggr(:,:,:,:) 3519 real(dp),allocatable :: symrec_cart(:,:,:) 3520 real(dp),pointer :: grad(:,:,:) 3521 type(coeff3_type),target,allocatable :: tmp_grhoij(:) 3522 type(pawrhoij_type),pointer :: pawrhoij_unsym_all(:) 3523 3524 ! ********************************************************************* 3525 3526 !Sizes of pawrhoij datastructures 3527 nrhoij=size(pawrhoij) 3528 nrhoij_unsym=size(pawrhoij_unsym) 3529 3530 !Set up parallelism over atoms 3531 paral_atom=(present(comm_atom).and.(nrhoij/=natom)) 3532 paral_atom_unsym=(present(comm_atom).and.(nrhoij_unsym/=natom)) 3533 nullify(my_atmtab);if (present(mpi_atmtab)) my_atmtab => mpi_atmtab 3534 my_comm_atom=xmpi_comm_self;if (present(comm_atom)) my_comm_atom=comm_atom 3535 call get_my_atmtab(my_comm_atom,my_atmtab,my_atmtab_allocated,paral_atom,natom) 3536 3537 !Test: consistency between choice and ngrhoij 3538 ngrhoij=0 3539 if (nrhoij>0) then 3540 ngrhoij=pawrhoij(1)%ngrhoij 3541 if(choice==2) ngrhoij=min(3,pawrhoij(1)%ngrhoij) 3542 if(choice==3.or.choice==4) ngrhoij=min(6,pawrhoij(1)%ngrhoij) 3543 if(choice==23.or.choice==24) ngrhoij=min(9,pawrhoij(1)%ngrhoij) 3544 if ((choice==1.and.ngrhoij/=0).or.(choice==2.and.ngrhoij/=3).or. & 3545 & (choice==3.and.ngrhoij/=6).or.(choice==23.and.ngrhoij/=9).or. & 3546 & (choice==4.and.ngrhoij/=6).or.(choice==24.and.ngrhoij/=9) ) then 3547 msg='Inconsistency between variables choice and ngrhoij !' 3548 LIBPAW_BUG(msg) 3549 end if 3550 end if 3551 3552 !Antiferro case ? 3553 antiferro=.false.;if (nrhoij>0) antiferro=(pawrhoij(1)%nspden==2.and.pawrhoij(1)%nsppol==1) 3554 !Non-collinear case 3555 noncoll=.false.;if (nrhoij>0) noncoll=(pawrhoij(1)%nspden==4) 3556 !Do we use antiferro symmetries ? 3557 use_afm=((antiferro).or.(noncoll.and.afm_noncoll)) 3558 !Do we impose zero magnetization? 3559 use_zeromag_=.false. ; if (present(use_zeromag)) use_zeromag_=use_zeromag 3560 3561 ! Does not symmetrize imaginary part for GS calculations 3562 cplex_eff=1 3563 if (nrhoij>0.and.(ipert>0.or.antiferro.or.noncoll)) cplex_eff=pawrhoij(1)%cplex_rhoij 3564 3565 !Do we have a phase due to q-vector? 3566 has_qphase=.false. 3567 if (nrhoij>0) then 3568 has_qphase=(pawrhoij(1)%qphase==2) 3569 if (present(qphon)) then 3570 if (any(abs(qphon(1:3))>tol8).and.(.not.has_qphase)) then 3571 msg='Should have qphase=2 for a non-zero q!' 3572 LIBPAW_BUG(msg) 3573 end if 3574 end if 3575 end if 3576 3577 !Printing of unsymetrized Rhoij 3578 if (nrhoij>0.and.optrhoij==1.and.pawprtvol/=-10001) then 3579 wrt_mode='COLL';if (paral_atom) wrt_mode='PERS' 3580 pertstrg="RHOIJ";if (ipert>0) pertstrg="RHOIJ(1)" 3581 natinc=1;if(nrhoij>1.and.pawprtvol>=0) natinc=nrhoij-1 3582 write(msg, '(7a)') ch10," PAW TEST:",ch10,& 3583 & ' ========= Values of ',trim(pertstrg),' before symetrization =========',ch10 3584 call wrtout(std_out,msg,wrt_mode) 3585 do iatm=1,nrhoij,natinc 3586 iatom=iatm; if (paral_atom) iatom=my_atmtab(iatm) 3587 if (nrhoij==1.and.ipert>0.and.ipert<=natom) iatom=ipert 3588 call pawrhoij_print_rhoij(pawrhoij_unsym(iatm)%rhoij_,pawrhoij_unsym(iatm)%cplex_rhoij,& 3589 & pawrhoij_unsym(iatm)%qphase,iatom,natom,& 3590 & unit=std_out,opt_prtvol=pawprtvol,mode_paral=wrt_mode) 3591 end do 3592 call wrtout(std_out,"",wrt_mode) 3593 end if 3594 3595 !Symetrization occurs only when nsym>1 3596 if (nsym>1) then 3597 3598 ! Symetrization of gradients not compatible with nspden=4 3599 if (nrhoij>0) then 3600 if (choice>2.and.pawrhoij(1)%nspden==4) then 3601 msg='For the time being, choice>2 is not compatible with nspden=4 !' 3602 LIBPAW_BUG(msg) 3603 end if 3604 end if 3605 3606 ! Symetry matrixes must be in memory 3607 if (pawang%nsym==0) then 3608 msg='pawang%zarot must be allocated !' 3609 LIBPAW_BUG(msg) 3610 end if 3611 3612 if (has_qphase.and.choice>1) then 3613 msg='choice>1 not compatible with q-phase !' 3614 LIBPAW_BUG(msg) 3615 end if 3616 3617 ! Several inits/allocations 3618 if (noncoll) then 3619 LIBPAW_ALLOCATE(symrec_cart,(3,3,nsym)) 3620 do irot=1,nsym 3621 symrec_cart(:,:,irot)=symrhoij_symcart(gprimd,rprimd,symrec(:,:,irot)) 3622 end do 3623 end if 3624 ishift2=0;ishift3=0;ishift4=0 3625 if (choice>1) then 3626 iafm=merge(2,1,antiferro) 3627 qphase=merge(2,1,has_qphase) 3628 LIBPAW_ALLOCATE(sumgr,(cplex_eff,ngrhoij,qphase)) 3629 LIBPAW_ALLOCATE(rotgr,(cplex_eff,ngrhoij,iafm,qphase)) 3630 if (noncoll) then 3631 LIBPAW_ALLOCATE(summaggr,(cplex_eff,ngrhoij,3,qphase)) 3632 LIBPAW_ALLOCATE(rotmaggr,(cplex_eff,ngrhoij,3,qphase)) 3633 end if 3634 if (choice==23) ishift2=6 3635 if (choice==24) ishift4=3 3636 if (.not.paral_atom_unsym) then 3637 ! Have to make a temporary copy of grhoij 3638 LIBPAW_DATATYPE_ALLOCATE(tmp_grhoij,(nrhoij)) 3639 do iatm=1,nrhoij 3640 sz1=pawrhoij_unsym(iatm)%cplex_rhoij*pawrhoij_unsym(iatm)%qphase & 3641 & *pawrhoij_unsym(iatm)%lmn2_size 3642 sz2=pawrhoij_unsym(iatm)%nspden 3643 LIBPAW_ALLOCATE(tmp_grhoij(iatm)%value,(ngrhoij,sz1,sz2)) 3644 tmp_grhoij(iatm)%value(1:ngrhoij,1:sz1,1:sz2)= & 3645 & pawrhoij_unsym(iatm)%grhoij(1:ngrhoij,1:sz1,1:sz2) 3646 end do 3647 end if 3648 end if 3649 3650 ! In case of paw_rhoij_unsym distributed over atomic sites, gather it 3651 if (paral_atom_unsym) then 3652 LIBPAW_DATATYPE_ALLOCATE(pawrhoij_unsym_all,(natom)) 3653 call pawrhoij_nullify(pawrhoij_unsym_all) 3654 call pawrhoij_gather(pawrhoij_unsym,pawrhoij_unsym_all,-1,my_comm_atom,& 3655 & with_lmnmix=.false.,with_rhoijp=.false.,& 3656 & with_rhoijres=.false.,with_grhoij=(choice>1)) 3657 nrhoij1=natom 3658 else 3659 pawrhoij_unsym_all=>pawrhoij_unsym 3660 nrhoij1=nrhoij_unsym 3661 end if 3662 3663 3664 ! Loops over atoms and spin components 3665 ! ------------------------------------ 3666 do iatm=1,nrhoij 3667 iatom=iatm;if (paral_atom) iatom=my_atmtab(iatm) 3668 if (nrhoij==1.and.ipert>0.and.ipert<=natom.and.(.not.paral_atom)) iatom=ipert 3669 itypat=typat(iatom) 3670 lmn_size=pawrhoij(iatm)%lmn_size 3671 lmn2_size=pawrhoij(iatm)%lmn2_size 3672 qphase=pawrhoij(iatm)%qphase 3673 cplex_rhoij=pawrhoij(iatm)%cplex_rhoij 3674 cplex_eff=1;if (ipert>0.or.antiferro.or.noncoll) cplex_eff=cplex_rhoij 3675 use_res=(pawrhoij(iatm)%use_rhoijres>0) 3676 indlmn => pawtab(itypat)%indlmn 3677 3678 nselect=0 3679 do ispden=1,pawrhoij(iatm)%nsppol 3680 jspden=min(3-ispden,pawrhoij(iatm)%nsppol) 3681 3682 ! Store old -rhoij in residual 3683 if (optrhoij==1.and.use_res) then 3684 pawrhoij(iatm)%rhoijres(:,ispden)=zero 3685 if (noncoll) pawrhoij(iatm)%rhoijres(:,2:4)=zero 3686 if (antiferro) pawrhoij(iatm)%rhoijres(:,2)=zero 3687 do iq=1,qphase 3688 iq0=(iq-1)*cplex_rhoij*lmn2_size 3689 if (cplex_rhoij==1) then 3690 do irhoij=1,pawrhoij(iatm)%nrhoijsel 3691 klmn1=iq0+pawrhoij(iatm)%rhoijselect(irhoij);jrhoij=iq0+irhoij 3692 pawrhoij(iatm)%rhoijres(klmn1,ispden)=-pawrhoij(iatm)%rhoijp(jrhoij,ispden) 3693 end do 3694 else 3695 do irhoij=1,pawrhoij(iatm)%nrhoijsel 3696 klmn1=iq0+2*pawrhoij(iatm)%rhoijselect(irhoij);jrhoij=iq0+2*irhoij 3697 pawrhoij(iatm)%rhoijres(klmn1-1:klmn1,ispden)=-pawrhoij(iatm)%rhoijp(jrhoij-1:jrhoij,ispden) 3698 end do 3699 end if 3700 if (noncoll) then 3701 if (cplex_rhoij==1) then 3702 do mu=2,4 3703 do irhoij=1,pawrhoij(iatm)%nrhoijsel 3704 klmn1=iq0+pawrhoij(iatm)%rhoijselect(irhoij);jrhoij=iq0+irhoij 3705 pawrhoij(iatm)%rhoijres(klmn1,mu)=-pawrhoij(iatm)%rhoijp(jrhoij,mu) 3706 end do 3707 end do 3708 else 3709 do mu=2,4 3710 do irhoij=1,pawrhoij(iatm)%nrhoijsel 3711 klmn1=iq0+2*pawrhoij(iatm)%rhoijselect(irhoij);jrhoij=iq0+2*irhoij 3712 pawrhoij(iatm)%rhoijres(klmn1-1:klmn1,mu)=-pawrhoij(iatm)%rhoijp(jrhoij-1:jrhoij,mu) 3713 end do 3714 end do 3715 end if 3716 end if 3717 if (antiferro) then 3718 if (cplex_rhoij==1) then 3719 do irhoij=1,pawrhoij(iatm)%nrhoijsel 3720 klmn1=iq0+pawrhoij(iatm)%rhoijselect(irhoij);jrhoij=iq0+irhoij 3721 pawrhoij(iatm)%rhoijres(klmn1,2)=-pawrhoij(iatm)%rhoijp(jrhoij,2) 3722 end do 3723 else 3724 do irhoij=1,pawrhoij(iatm)%nrhoijsel 3725 klmn1=iq0+2*pawrhoij(iatm)%rhoijselect(irhoij);jrhoij=iq0+2*irhoij 3726 pawrhoij(iatm)%rhoijres(klmn1-1:klmn1,2)=-pawrhoij(iatm)%rhoijp(jrhoij-1:jrhoij,2) 3727 end do 3728 end if 3729 end if 3730 end do 3731 end if 3732 3733 3734 ! Loops over (il,im) and (jl,jm) 3735 ! ------------------------------ 3736 jl0=-1;jln0=-1;indexj=1 3737 do jlmn=1,lmn_size 3738 jl=indlmn(1,jlmn) 3739 jlpm=1+jl+indlmn(2,jlmn) 3740 jln=indlmn(5,jlmn) 3741 if (jln/=jln0) indexj=indexj+2*jl0+1 3742 j0lmn=jlmn*(jlmn-1)/2 3743 il0=-1;iln0=-1;indexi=1 3744 do ilmn=1,jlmn 3745 il=indlmn(1,ilmn) 3746 ilpm=1+il+indlmn(2,ilmn) 3747 iln=indlmn(5,ilmn) 3748 if (iln/=iln0) indexi=indexi+2*il0+1 3749 klmn=j0lmn+ilmn 3750 klmn1=merge(klmn,2*klmn-1,cplex_rhoij==1) 3751 3752 nsym_used(:)=0 3753 if (optrhoij==1) rotrho=zero 3754 if (optrhoij==1.and.noncoll) rotmag=zero 3755 if (choice>1) rotgr=zero 3756 if (choice>1.and.noncoll) rotmaggr=zero 3757 3758 3759 ! Loop over symmetries 3760 ! -------------------- 3761 do irot=1,nsym 3762 3763 if ((symafm(irot)/=1).and.(.not.use_afm)) cycle 3764 kspden=ispden;if (symafm(irot)==-1) kspden=jspden 3765 iafm=1;if ((antiferro).and.(symafm(irot)==-1)) iafm=2 3766 factafm=dble(symafm(irot)) 3767 3768 nsym_used(iafm)=nsym_used(iafm)+1 3769 at_indx=min(indsym(4,irot,iatom),nrhoij1) 3770 3771 if (has_qphase) then 3772 arg=two_pi*(qphon(1)*indsym(1,irot,iatom)+qphon(2)*indsym(2,irot,iatom) & 3773 & +qphon(3)*indsym(3,irot,iatom)) 3774 phase(1)=cos(arg);phase(2)=sin(arg) 3775 end if 3776 3777 if (optrhoij==1) sumrho=zero 3778 if (optrhoij==1.and.noncoll) summag=zero 3779 if (choice>1) sumgr=zero 3780 if (choice>1.and.noncoll) summaggr=zero 3781 3782 if (choice>1) then 3783 if (paral_atom_unsym) then 3784 grad => pawrhoij_unsym_all(at_indx)%grhoij 3785 else 3786 grad => tmp_grhoij(at_indx)%value 3787 end if 3788 end if 3789 3790 3791 ! Accumulate values over (mi,mj) 3792 ! ------------------------------ 3793 do mj=1,2*jl+1 3794 indexjj=indexj+mj;indexjj0=indexjj*(indexjj-1)/2 3795 do mi=1,2*il+1 3796 factsym(:)=one 3797 indexii=indexi+mi 3798 if (indexii<=indexjj) then 3799 indexk=indexjj0+indexii 3800 factsym(2)=one 3801 else 3802 indexk=indexii*(indexii-1)/2+indexjj 3803 factsym(2)=-one 3804 end if 3805 indexkc=cplex_rhoij*(indexk-1) 3806 indexkc_q=indexkc+cplex_rhoij*lmn2_size 3807 3808 ! Be careful: use here R_rel^-1 in term of spherical harmonics 3809 ! which is tR_rec in term of spherical harmonics 3810 ! so, use transpose[zarot] 3811 zarot2=pawang%zarot(mi,ilpm,il+1,irot)*pawang%zarot(mj,jlpm,jl+1,irot) 3812 ! zarot2=pawang%zarot(ilpm,mi,il+1,irot)*pawang%zarot(jlpm,mj,jl+1,irot) 3813 3814 ! Rotate rhoij 3815 if (optrhoij==1) then 3816 fact(1)=factsym(1);fact(2)=factsym(2)*factafm !????? What? MT 3817 sumrho(1:cplex_eff,iafm,1)=sumrho(1:cplex_eff,iafm,1) & 3818 & +fact(1:cplex_eff)*zarot2 & 3819 & *pawrhoij_unsym_all(at_indx)%rhoij_(indexkc+1:indexkc+cplex_eff,kspden) 3820 if (qphase==2) & 3821 & sumrho(1:cplex_eff,iafm,2)=sumrho(1:cplex_eff,iafm,2) & 3822 & +fact(1:cplex_eff)*zarot2 & 3823 & *pawrhoij_unsym_all(at_indx)%rhoij_(indexkc_q+1:indexkc_q+cplex_eff,kspden) 3824 end if 3825 3826 ! Rotate rhoij magnetization 3827 if (optrhoij==1.and.noncoll) then 3828 fact(1)=factsym(1)*factafm;fact(2)=factsym(2) 3829 do mu=1,3 3830 summag(1:cplex_eff,mu,1)=summag(1:cplex_eff,mu,1) & 3831 & +fact(1:cplex_eff)*zarot2 & 3832 & *pawrhoij_unsym_all(at_indx)%rhoij_(indexkc+1:indexkc+cplex_eff,1+mu) 3833 end do 3834 if (qphase==2) then 3835 do mu=1,3 3836 summag(1:cplex_eff,mu,2)=summag(1:cplex_eff,mu,2) & 3837 & +fact(1:cplex_eff)*zarot2 & 3838 & *pawrhoij_unsym_all(at_indx)%rhoij_(indexkc_q+1:indexkc_q+cplex_eff,1+mu) 3839 end do 3840 end if 3841 end if 3842 3843 ! Rotate gradients of rhoij 3844 if (choice>1) then 3845 fact(1)=factsym(1);fact(2)=factsym(2)*factafm !????? What? MT 3846 do mu=1,ngrhoij 3847 sumgr(1:cplex_eff,mu,1)=sumgr(1:cplex_eff,mu,1) & 3848 & +fact(1:cplex_eff)*zarot2*grad(mu,indexkc+1:indexkc+cplex_eff,kspden) 3849 end do 3850 if (qphase==2) then 3851 do mu=1,ngrhoij 3852 sumgr(1:cplex_eff,mu,2)=sumgr(1:cplex_eff,mu,2) & 3853 & +fact(1:cplex_eff)*zarot2*grad(mu,indexkc_q+1:indexkc_q+cplex_eff,kspden) 3854 end do 3855 end if 3856 end if 3857 3858 ! Rotate gradients of rhoij magnetization 3859 if (choice>1.and.noncoll) then 3860 fact(1)=factsym(1)*factafm;fact(2)=factsym(2) 3861 do mu=1,3 3862 do nu=1,ngrhoij 3863 summaggr(1:cplex_eff,nu,mu,1)=summaggr(1:cplex_eff,nu,mu,1) & 3864 & +fact(1:cplex_eff)*zarot2*grad(nu,indexkc+1:indexkc+cplex_eff,1+mu) 3865 end do 3866 end do 3867 if (qphase==2) then 3868 do mu=1,3 3869 do nu=1,ngrhoij 3870 summaggr(1:cplex_eff,nu,mu,2)=summaggr(1:cplex_eff,nu,mu,2) & 3871 & +fact(1:cplex_eff)*zarot2*grad(nu,indexkc_q+1:indexkc_q+cplex_eff,1+mu) 3872 end do 3873 end do 3874 end if 3875 end if 3876 3877 end do ! mi 3878 end do ! mj 3879 3880 ! Apply phase for phonons 3881 if (has_qphase) then 3882 !Remember, RHOij is stored as follows: 3883 ! RHOij= [rhoij(2klmn-1)+i.rhoij(2klmn)] 3884 ! +i.[rhoij(lnm2_size+2klmn-1)+i.rhoij(lmn2_size+2klmn)] 3885 if (optrhoij==1) then 3886 do iplex=1,cplex_rhoij 3887 rhoijc(1)=sumrho(iplex,iafm,1) 3888 rhoijc(2)=sumrho(iplex,iafm,2) 3889 sumrho(iplex,iafm,1)=phase(1)*rhoijc(1)-phase(2)*rhoijc(2) 3890 sumrho(iplex,iafm,2)=phase(1)*rhoijc(2)+phase(2)*rhoijc(1) 3891 end do 3892 if (noncoll) then 3893 do iplex=1,cplex_rhoij 3894 do mu=1,3 3895 rhoijc(1)=summag(iplex,mu,1) 3896 rhoijc(2)=summag(iplex,mu,2) 3897 summag(iplex,mu,1)=phase(1)*rhoijc(1)-phase(2)*rhoijc(2) 3898 summag(iplex,mu,2)=phase(1)*rhoijc(2)+phase(2)*rhoijc(1) 3899 end do 3900 end do 3901 end if 3902 end if 3903 if (choice>1) then 3904 do iplex=1,cplex_rhoij 3905 do mu=1,3 3906 rhoijc(1)=sumgr(iplex,mu,1) 3907 rhoijc(2)=sumgr(iplex,mu,2) 3908 sumgr(iplex,mu,1)=phase(1)*rhoijc(1)-phase(2)*rhoijc(2) 3909 sumgr(iplex,mu,2)=phase(1)*rhoijc(2)+phase(2)*rhoijc(1) 3910 end do 3911 end do 3912 if (noncoll) then 3913 do iplex=1,cplex_rhoij 3914 do mu=1,3 3915 do nu=1,ngrhoij 3916 rhoijc(1)=summaggr(iplex,nu,mu,1) 3917 rhoijc(2)=summaggr(iplex,nu,mu,2) 3918 summaggr(iplex,nu,mu,1)=phase(1)*rhoijc(1)-phase(2)*rhoijc(2) 3919 summaggr(iplex,nu,mu,2)=phase(1)*rhoijc(2)+phase(2)*rhoijc(1) 3920 end do 3921 end do 3922 end do 3923 end if 3924 end if 3925 end if 3926 3927 ! Add contribution of this rotation 3928 if (optrhoij==1) then 3929 do iq=1,qphase 3930 rotrho(1:cplex_eff,iafm,iq)=rotrho(1:cplex_eff,iafm,iq) & 3931 & +sumrho(1:cplex_eff,iafm,iq) 3932 end do 3933 end if 3934 3935 3936 ! Rotate vector fields in real space (forces, magnetization, etc...) 3937 ! Should use symrel^-1 but use transpose[symrec] instead 3938 ! ===== Rhoij magnetization ==== 3939 if (noncoll.and.optrhoij==1) then 3940 do iq=1,qphase 3941 do nu=1,3 3942 do mu=1,3 3943 rotmag(1:cplex_eff,mu,iq)=rotmag(1:cplex_eff,mu,iq) & 3944 & +symrec_cart(mu,nu,irot)*summag(1:cplex_eff,nu,iq) 3945 end do 3946 end do 3947 end do 3948 end if 3949 ! ===== Derivatives vs atomic positions ==== 3950 if (choice==2.or.choice==23.or.choice==24) then 3951 do iq=1,qphase 3952 do nu=1,3 3953 nushift=nu+ishift2 3954 do mu=1,3 3955 mushift=mu+ishift2 3956 rotgr(1:cplex_eff,mushift,iafm,iq)=rotgr(1:cplex_eff,mushift,iafm,iq) & 3957 & +dble(symrec(mu,nu,irot))*sumgr(1:cplex_eff,nushift,iq) 3958 end do 3959 end do 3960 end do 3961 if (noncoll) then 3962 do iq=1,qphase 3963 do mub=1,3 ! Loop on magnetization components 3964 do mua=1,3 ! Loop on gradients 3965 mushift=mua+ishift2 3966 sum1(:)=zero;xsym(1:3)=dble(symrec(mua,1:3,irot)) 3967 do nu=1,3 3968 syma=symrec_cart(mub,nu,irot) 3969 sum1(1:cplex_eff)=sum1(1:cplex_eff)+syma & 3970 & *(summaggr(1:cplex_eff,ishift2+1,nu,iq)*xsym(1) & 3971 & +summaggr(1:cplex_eff,ishift2+2,nu,iq)*xsym(2) & 3972 & +summaggr(1:cplex_eff,ishift2+3,nu,iq)*xsym(3)) 3973 end do 3974 rotmaggr(1:cplex_eff,mushift,mub,iq)= & 3975 & rotmaggr(1:cplex_eff,mushift,mub,iq)+sum1(1:cplex_eff) 3976 end do 3977 end do 3978 end do 3979 end if 3980 end if 3981 ! ===== Derivatives vs strain ==== 3982 if (choice==3.or.choice==23) then 3983 do iq=1,qphase 3984 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) 3985 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) 3986 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) 3987 work1(1:cplex_eff,3,1)=work1(1:cplex_eff,1,3);work1(1:cplex_eff,3,2)=work1(1:cplex_eff,2,3) 3988 work1(1:cplex_eff,2,1)=work1(1:cplex_eff,1,2) 3989 do mu=1,6 3990 mushift=mu+ishift3 3991 mua=alpha(mu);mub=beta(mu) 3992 sum1(:)=zero;xsym(1:3)=dble(symrec(mub,1:3,irot)) 3993 do nu=1,3 3994 syma=dble(symrec(mua,nu,irot)) 3995 sum1(1:cplex_eff)=sum1(1:cplex_eff) & 3996 & +syma*(work1(1:cplex_eff,nu,1)*xsym(1) & 3997 & +work1(1:cplex_eff,nu,2)*xsym(2) & 3998 & +work1(1:cplex_eff,nu,3)*xsym(3)) 3999 end do 4000 rotgr(1:cplex_eff,mushift,iafm,iq)= & 4001 & rotgr(1:cplex_eff,mushift,iafm,iq)+sum1(1:cplex_eff) 4002 end do 4003 end do 4004 end if 4005 ! ===== Second derivatives vs atomic positions ==== 4006 if (choice==4.or.choice==24) then 4007 do iq=1,qphase 4008 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) 4009 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) 4010 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) 4011 work1(1:cplex_eff,3,1)=work1(1:cplex_eff,1,3);work1(1:cplex_eff,3,2)=work1(1:cplex_eff,2,3) 4012 work1(1:cplex_eff,2,1)=work1(1:cplex_eff,1,2) 4013 do mu=1,6 4014 mushift=mu+ishift4 4015 mua=alpha(mu);mub=beta(mu) 4016 sum1(:)=zero 4017 xsym(1:3)=dble(symrec(mub,1:3,irot)) 4018 do nu=1,3 4019 syma=dble(symrec(mua,nu,irot)) 4020 sum1(1:cplex_eff)=sum1(1:cplex_eff) & 4021 & +syma*(work1(1:cplex_eff,nu,1)*xsym(1) & 4022 & +work1(1:cplex_eff,nu,2)*xsym(2) & 4023 & +work1(1:cplex_eff,nu,3)*xsym(3)) 4024 end do 4025 rotgr(1:cplex_eff,mushift,iafm,iq)= & 4026 & rotgr(1:cplex_eff,mushift,iafm,iq)+sum1(1:cplex_eff) 4027 end do 4028 end do 4029 end if 4030 4031 end do ! End loop over symmetries 4032 4033 4034 ! Store average result (over symmetries) 4035 ! -------------------------------------- 4036 4037 ! Rhoij 4038 if (optrhoij==1) then 4039 do iq=1,qphase 4040 klmn1q=klmn1+(iq-1)*lmn2_size 4041 pawrhoij(iatm)%rhoijp(klmn1q,ispden)=rotrho(1,1,iq)/nsym_used(1) 4042 if (cplex_rhoij==2) then 4043 if (cplex_eff==1) ro=pawrhoij_unsym_all(iatom)%rhoij_(klmn1q+1,ispden) 4044 if (cplex_eff==2) ro=rotrho(2,1,iq)/nsym_used(1) 4045 pawrhoij(iatm)%rhoijp(klmn1q+1,ispden)=ro 4046 end if 4047 end do 4048 end if 4049 4050 ! Rhoij magnetization 4051 if (noncoll.and.optrhoij==1) then 4052 do mu=2,4 4053 do iq=1,qphase 4054 klmn1q=klmn1+(iq-1)*lmn2_size 4055 if (use_zeromag_) then 4056 pawrhoij(iatm)%rhoijp(klmn1q,mu)=zero 4057 else 4058 pawrhoij(iatm)%rhoijp(klmn1q,mu)=rotmag(1,mu-1,iq)/nsym_used(1) 4059 end if 4060 if (cplex_rhoij==2) then 4061 if (cplex_eff==1) ro=pawrhoij_unsym_all(iatom)%rhoij_(klmn1q+1,mu) 4062 if (cplex_eff==2) ro=rotmag(2,mu-1,iq)/nsym_used(1) 4063 pawrhoij(iatm)%rhoijp(klmn1q+1,mu)=ro 4064 end if 4065 end do 4066 end do 4067 end if 4068 4069 ! Rhoij^down when antiferro 4070 if (antiferro.and.optrhoij==1) then 4071 if (nsym_used(2)>0) then 4072 do iq=1,qphase 4073 klmn1q=klmn1+(iq-1)*lmn2_size 4074 pawrhoij(iatm)%rhoijp(klmn1q,2)=rotrho(1,2,iq)/nsym_used(2) 4075 if (cplex_rhoij==2) then 4076 if (cplex_eff==1) ro=pawrhoij_unsym_all(iatom)%rhoij_(klmn1q+1,2) 4077 if (cplex_eff==2) ro=rotrho(2,2,iq)/nsym_used(2) 4078 pawrhoij(iatm)%rhoijp(klmn1q+1,2)=ro 4079 end if 4080 end do 4081 end if 4082 end if 4083 4084 ! Gradients of rhoij 4085 if (choice>1) then 4086 do iq=1,qphase 4087 klmn1q=klmn1+(iq-1)*lmn2_size 4088 do iplex=1,cplex_eff 4089 do mu=1,ngrhoij 4090 pawrhoij(iatm)%grhoij(mu,klmn1q,ispden)=rotgr(iplex,mu,1,iq)/nsym_used(1) 4091 end do 4092 if (noncoll) then 4093 if (use_zeromag_.and.iplex==1) then 4094 pawrhoij(iatm)%grhoij(mu,klmn1q,2:4)=zero 4095 else 4096 do nu=1,3 4097 pawrhoij(iatm)%grhoij(mu,klmn1q,1+nu)=rotmaggr(iplex,mu,nu,iq)/nsym_used(1) 4098 end do 4099 end if 4100 end if 4101 if (antiferro.and.nsym_used(2)>0) then 4102 do mu=1,ngrhoij 4103 pawrhoij(iatm)%grhoij(mu,klmn1q,ispden)=rotgr(iplex,mu,2,iq)/nsym_used(2) 4104 end do 4105 end if 4106 klmn1q=klmn1q+1 4107 end do 4108 !if cplex_eff<cplex_rhoij, imaginary part of grhoij is unchanged 4109 end do 4110 end if 4111 4112 4113 il0=il;iln0=iln ! End loops over (il,im) and (jl,jm) 4114 end do 4115 jl0=jl;jln0=jln 4116 end do 4117 4118 end do ! End loop over ispden 4119 4120 ! Select non-zero elements of rhoij 4121 if (optrhoij==1) then 4122 call pawrhoij_filter(pawrhoij(iatm)%rhoijp,pawrhoij(iatm)%rhoijselect,& 4123 & pawrhoij(iatm)%nrhoijsel,cplex_rhoij,qphase,lmn2_size,& 4124 & pawrhoij(iatm)%nspden) 4125 end if 4126 4127 ! Add new rhoij to rhoij residual 4128 if (optrhoij==1.and.use_res) then 4129 do ispden=1,pawrhoij(iatm)%nspden 4130 do iq=1,qphase 4131 iq0=(iq-1)*lmn2_size 4132 if (cplex_rhoij==1) then 4133 do irhoij=1,pawrhoij(iatm)%nrhoijsel 4134 klmn1=iq0+pawrhoij(iatm)%rhoijselect(irhoij) ; jrhoij=iq0+irhoij 4135 pawrhoij(iatm)%rhoijres(klmn1,ispden)= & 4136 & pawrhoij(iatm)%rhoijres(klmn1,ispden) & 4137 & +pawrhoij(iatm)%rhoijp(jrhoij,ispden) 4138 end do 4139 else 4140 do irhoij=1,pawrhoij(iatm)%nrhoijsel 4141 klmn1=iq0+2*pawrhoij(iatm)%rhoijselect(irhoij) ; jrhoij=iq0+2*irhoij 4142 pawrhoij(iatm)%rhoijres(klmn1-1:klmn1,ispden)= & 4143 & pawrhoij(iatm)%rhoijres(klmn1-1:klmn1,ispden) & 4144 & +pawrhoij(iatm)%rhoijp(jrhoij-1:jrhoij,ispden) 4145 end do 4146 end if 4147 end do 4148 end do 4149 end if 4150 4151 end do ! End loop over atoms 4152 4153 if (noncoll) then 4154 LIBPAW_DEALLOCATE(symrec_cart) 4155 end if 4156 if (choice>1) then 4157 if (.not.paral_atom_unsym) then 4158 do iatm=1,nrhoij 4159 LIBPAW_DEALLOCATE(tmp_grhoij(iatm)%value) 4160 end do 4161 LIBPAW_DATATYPE_DEALLOCATE(tmp_grhoij) 4162 end if 4163 LIBPAW_DEALLOCATE(sumgr) 4164 LIBPAW_DEALLOCATE(rotgr) 4165 if (noncoll) then 4166 LIBPAW_DEALLOCATE(summaggr) 4167 LIBPAW_DEALLOCATE(rotmaggr) 4168 end if 4169 end if 4170 if(paral_atom_unsym) then 4171 call pawrhoij_free(pawrhoij_unsym_all) 4172 LIBPAW_DATATYPE_DEALLOCATE(pawrhoij_unsym_all) 4173 end if 4174 4175 4176 else ! nsym>1 4177 4178 ! ********************************************************************* 4179 ! If nsym==1, only copy rhoij_ into rhoij 4180 ! also has to fill rhoijselect array 4181 4182 if (antiferro) then 4183 msg=' In the antiferromagnetic case, nsym cannot be 1' 4184 LIBPAW_BUG(msg) 4185 end if 4186 4187 if (optrhoij==1) then 4188 4189 do iatm=1,nrhoij 4190 iatom=iatm;if ((paral_atom).and.(.not.paral_atom_unsym)) iatom=my_atmtab(iatm) 4191 cplex_rhoij=pawrhoij(iatm)%cplex_rhoij 4192 qphase=pawrhoij(iatm)%qphase 4193 lmn2_size=pawrhoij(iatm)%lmn2_size 4194 use_res=(pawrhoij(iatm)%use_rhoijres>0) 4195 4196 ! Store -rhoij_input in rhoij residual 4197 if (use_res) then 4198 pawrhoij(iatm)%rhoijres(:,:)=zero 4199 do iq=1,qphase 4200 iq0=(iq-1)*lmn2_size 4201 if (cplex_rhoij==1) then 4202 do ispden=1,pawrhoij(iatm)%nspden 4203 do irhoij=1,pawrhoij(iatm)%nrhoijsel 4204 klmn=iq0+pawrhoij(iatm)%rhoijselect(irhoij);jrhoij=iq0+irhoij 4205 pawrhoij(iatm)%rhoijres(klmn,ispden)=-pawrhoij(iatm)%rhoijp(jrhoij,ispden) 4206 end do 4207 end do 4208 else 4209 do ispden=1,pawrhoij(iatm)%nspden 4210 do irhoij=1,pawrhoij(iatm)%nrhoijsel 4211 klmn1=iq0+2*pawrhoij(iatm)%rhoijselect(irhoij);jrhoij=iq0+2*irhoij 4212 pawrhoij(iatm)%rhoijres(klmn1-1,ispden)=-pawrhoij(iatm)%rhoijp(jrhoij-1,ispden) 4213 pawrhoij(iatm)%rhoijres(klmn1 ,ispden)=-pawrhoij(iatm)%rhoijp(jrhoij ,ispden) 4214 end do 4215 end do 4216 end if 4217 end do 4218 end if 4219 4220 ! Select non-zero elements of input rhoij 4221 call pawrhoij_filter(pawrhoij(iatm)%rhoijp,pawrhoij(iatm)%rhoijselect,& 4222 & pawrhoij(iatm)%nrhoijsel,cplex_rhoij,qphase,lmn2_size,& 4223 & pawrhoij(iatm)%nspden,rhoij_input=pawrhoij_unsym(iatom)%rhoij_) 4224 4225 ! Add new rhoij to rhoij residual 4226 if (use_res) then 4227 do ispden=1,pawrhoij(iatm)%nspden 4228 do iq=1,qphase 4229 iq0=(iq-1)*lmn2_size 4230 if (cplex_rhoij==1) then 4231 do irhoij=1,pawrhoij(iatm)%nrhoijsel 4232 klmn1=iq0+pawrhoij(iatm)%rhoijselect(irhoij) ; jrhoij=iq0+irhoij 4233 pawrhoij(iatm)%rhoijres(klmn1,ispden)= & 4234 & pawrhoij(iatm)%rhoijres(klmn1,ispden) & 4235 & +pawrhoij(iatm)%rhoijp(jrhoij,ispden) 4236 end do 4237 else 4238 do irhoij=1,pawrhoij(iatm)%nrhoijsel 4239 klmn1=iq0+2*pawrhoij(iatm)%rhoijselect(irhoij) ; jrhoij=iq0+2*irhoij 4240 pawrhoij(iatm)%rhoijres(klmn1-1:klmn1,ispden)= & 4241 & pawrhoij(iatm)%rhoijres(klmn1-1:klmn1,ispden) & 4242 & +pawrhoij(iatm)%rhoijp(jrhoij-1:jrhoij,ispden) 4243 end do 4244 end if 4245 end do 4246 end do 4247 end if 4248 4249 end do ! iatm 4250 end if ! optrhoij 4251 4252 end if 4253 4254 !********************************************************************* 4255 !Printing of symetrized Rhoij 4256 if (nrhoij>0.and.optrhoij==1.and.pawprtvol/=-10001) then 4257 wrt_mode='COLL';if (paral_atom) wrt_mode='PERS' 4258 pertstrg="RHOIJ";if (ipert>0) pertstrg="RHOIJ(1)" 4259 natinc=1;if(nrhoij>1.and.pawprtvol>=0) natinc=nrhoij-1 4260 write(msg, '(7a)') ch10," PAW TEST:",ch10,& 4261 & ' ========= Values of ',trim(pertstrg),' after symetrization =========',ch10 4262 call wrtout(std_out,msg,wrt_mode) 4263 do iatm=1,nrhoij,natinc 4264 iatom=iatm; if (paral_atom) iatom=my_atmtab(iatm) 4265 if (nrhoij==1.and.ipert>0.and.ipert<=natom) iatom=ipert 4266 call pawrhoij_print_rhoij(pawrhoij(iatm)%rhoijp,pawrhoij(iatm)%cplex_rhoij,& 4267 & pawrhoij(iatm)%qphase,iatom,natom,& 4268 & rhoijselect=pawrhoij(iatm)%rhoijselect,unit=std_out,& 4269 & opt_prtvol=pawprtvol,mode_paral=wrt_mode) 4270 end do 4271 call wrtout(std_out,"",wrt_mode) 4272 end if 4273 4274 !Destroy atom table used for parallelism 4275 call free_my_atmtab(my_atmtab,my_atmtab_allocated) 4276 4277 !********************************************************************* 4278 !Small function: convert a symmetry operation 4279 !from reduced coordinates (integers) to cartesian coordinates (reals) 4280 contains 4281 function symrhoij_symcart(aprim,bprim,symred) 4282 4283 real(dp) :: symrhoij_symcart(3,3) 4284 integer,intent(in) :: symred(3,3) 4285 real(dp),intent(in) :: aprim(3,3),bprim(3,3) 4286 integer :: ii,jj,kk 4287 real(dp) :: tmp(3,3) 4288 symrhoij_symcart=zero;tmp=zero 4289 do kk=1,3 4290 do jj=1,3 4291 do ii=1,3 4292 tmp(ii,jj)=tmp(ii,jj)+bprim(ii,kk)*dble(symred(jj,kk)) 4293 end do 4294 end do 4295 end do 4296 do kk=1,3 4297 do jj=1,3 4298 do ii=1,3 4299 symrhoij_symcart(ii,jj)=symrhoij_symcart(ii,jj)+aprim(ii,kk)*tmp(jj,kk) 4300 end do 4301 end do 4302 end do 4303 end function symrhoij_symcart 4304 4305 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
2718 subroutine pawrhoij_unpack(rhoij) 2719 2720 !Arguments ------------------------------------ 2721 !scalars 2722 !arrays 2723 type(pawrhoij_type),intent(inout) :: rhoij(:) 2724 2725 !Local variables------------------------------- 2726 integer :: cplex_rhoij,natom,lmn2_size,nsp2,qphase 2727 integer :: i0,iat,iphase,isel,isppol,klmn 2728 2729 ! ************************************************************************* 2730 2731 natom = SIZE(rhoij) ; if (natom==0) return 2732 2733 do iat=1,natom 2734 2735 lmn2_size =rhoij(iat)%lmn2_size 2736 cplex_rhoij = rhoij(iat)%cplex_rhoij 2737 qphase = rhoij(iat)%qphase 2738 nsp2 = rhoij(iat)%nsppol;if (rhoij(iat)%nspden==4) nsp2=4 2739 2740 if (rhoij(iat)%use_rhoij_/=1) then ! Have to allocate rhoij 2741 LIBPAW_ALLOCATE(rhoij(iat)%rhoij_,(cplex_rhoij*qphase*lmn2_size,nsp2)) 2742 rhoij(iat)%use_rhoij_=1 2743 end if 2744 rhoij(iat)%rhoij_ = zero 2745 2746 do isppol=1,nsp2 2747 do iphase=1,qphase 2748 i0=(iphase-1)*lmn2_size*cplex_rhoij 2749 if (cplex_rhoij==1) then 2750 do isel=1,rhoij(iat)%nrhoijsel ! Looping over non-zero ij elements. 2751 klmn = rhoij(iat)%rhoijselect(isel) 2752 rhoij(iat)%rhoij_(i0+klmn,isppol) = rhoij(iat)%rhoijp(i0+isel,isppol) 2753 end do 2754 else 2755 do isel=1,rhoij(iat)%nrhoijsel ! Looping over non-zero ij elements. 2756 klmn = rhoij(iat)%rhoijselect(isel) 2757 rhoij(iat)%rhoij_(i0+2*klmn-1:i0+2*klmn,isppol) = & 2758 & rhoij(iat)%rhoijp(i0+2*isel-1:i0+2*isel,isppol) 2759 end do 2760 end if 2761 end do 2762 end do 2763 2764 end do ! natom 2765 2766 end subroutine pawrhoij_unpack