TABLE OF CONTENTS


m_paw_occupancies/initrhoij [ Functions ]

[ Top ] [ m_paw_occupancies ] [ Functions ]

NAME

 initrhoij

FUNCTION

 Initialize PAW rhoij occupancies (in packed storage)
 from atomic ones

INPUTS

  cplex=1 if rhoij are REAL, 2 if they are complex
  lexexch(ntypat)=l on which local exact-exchange is applied for a given type of atom
  lpawu(ntypat)=l on which U is applied for a given type of atom (PAW+U)
  mpi_atmtab(:)=--optional-- indexes of the atoms treated by current proc
  comm_atom=--optional-- MPI communicator over atoms
  my_natom=number of atoms treated by current processor
  natom=number of atoms
  nspden=number of spin-density components
  nspinor=number of spinorial components
  nsppol=1 for unpolarized, 2 for spin-polarized
  ntypat=number of atom types
  pawspnorb=flag: 1 if spin-orbit coupling is activated in PAW augmentation regions
  pawtab(ntypat) <type(pawtab_type)>=paw tabulated starting data
                                     (containing initial rhoij)
  spinat(3,natom)=initial spin of each atom, in unit of hbar/2.
  typat(natom)=type of each atom
  === Optional arguments
    ngrhoij=number of gradients to be allocated (OPTIONAL, default=0)
    nlmnmix=number of rhoij elements to be mixed during SCF cycle (OPTIONAL, default=0)
    use_rhoij_=1 if pawrhoij(:)%rhoij_ has to be allocated (OPTIONAL, default=0)
    use_rhoijres=1 if pawrhoij(:)%rhoijres has to be allocated (OPTIONAL, default=0)

OUTPUT

  pawrhoij(natom) <type(pawrhoij_type)>=rhoij quantities for each atom
                                        in packed storage

PARENTS

      gstate,respfn,setup_positron

CHILDREN

      free_my_atmtab,get_my_atmtab,pawrhoij_alloc

SOURCE

1066 subroutine initrhoij(cplex,lexexch,lpawu,my_natom,natom,&
1067 &                    nspden,nspinor,nsppol,ntypat,pawrhoij,pawspnorb,pawtab,spinat,typat,&
1068 &                    ngrhoij,nlmnmix,use_rhoij_,use_rhoijres,& ! optional arguments
1069 &                    mpi_atmtab,comm_atom) ! optional arguments (parallelism)
1070 
1071 
1072 !This section has been created automatically by the script Abilint (TD).
1073 !Do not modify the following lines by hand.
1074 #undef ABI_FUNC
1075 #define ABI_FUNC 'initrhoij'
1076 !End of the abilint section
1077 
1078  implicit none
1079 
1080 !Arguments ---------------------------------------------
1081 !scalars
1082  integer,intent(in) :: cplex,my_natom,natom,nspden,nspinor,nsppol,ntypat,pawspnorb
1083  integer,intent(in),optional :: comm_atom,ngrhoij,nlmnmix,use_rhoij_,use_rhoijres
1084  character(len=500) :: message
1085 !arrays
1086  integer,intent(in) :: lexexch(ntypat),lpawu(ntypat)
1087  integer,intent(in) :: typat(natom)
1088  integer,optional,target,intent(in) :: mpi_atmtab(:)
1089  real(dp),intent(in) :: spinat(3,natom)
1090  type(pawrhoij_type),intent(inout) :: pawrhoij(my_natom)
1091  type(pawtab_type),intent(in) :: pawtab(ntypat)
1092 
1093 !Local variables ---------------------------------------
1094 !Arrays
1095 !scalars
1096  integer :: iatom,iatom_rhoij,ilmn,ispden,itypat,j0lmn,jl,jlmn,jspden,klmn,klmn1,ln,lnspinat0,my_comm_atom
1097  integer :: ngrhoij0,nlmnmix0,nselect,nselect1,nspden_rhoij,use_rhoij_0,use_rhoijres0
1098  real(dp) :: ratio,ro,roshift,zratio,zz
1099  logical :: my_atmtab_allocated,paral_atom,spinat_zero,test_exexch,test_pawu,test_lnspinat
1100 !arrays
1101  integer,pointer :: my_atmtab(:),lnspinat(:)
1102  real(dp),allocatable :: occ(:)
1103 !************************************************************************
1104 
1105  DBG_ENTER("COLL")
1106 
1107 !PAW+U and local exact-exchange restriction
1108  do itypat=1,ntypat
1109    if (lpawu(itypat)/=lexexch(itypat).and. lpawu(itypat)/=-1.and.lexexch(itypat)/=-1) then
1110      message = ' lpawu must be equal to lexexch !'
1111      MSG_ERROR(message)
1112    end if
1113  end do
1114 
1115 !Set up parallelism over atoms
1116  paral_atom=(present(comm_atom).and.(my_natom/=natom))
1117  nullify(my_atmtab);if (present(mpi_atmtab)) my_atmtab => mpi_atmtab
1118  my_comm_atom=xmpi_comm_self;if (present(comm_atom)) my_comm_atom=comm_atom
1119  call get_my_atmtab(my_comm_atom,my_atmtab,my_atmtab_allocated,paral_atom,natom,my_natom_ref=my_natom)
1120 
1121  nspden_rhoij=pawrhoij_get_nspden(nspden,nspinor,pawspnorb)
1122  ratio=one;if (nspden_rhoij==2) ratio=half
1123  spinat_zero=all(abs(spinat(:,:))<tol10)
1124 
1125  if (my_natom>0) then
1126    ngrhoij0=0;if (present(ngrhoij)) ngrhoij0=ngrhoij
1127    nlmnmix0=0;if (present(nlmnmix)) nlmnmix0=nlmnmix
1128    use_rhoij_0=0;if (present(use_rhoij_)) use_rhoij_0=use_rhoij_
1129    use_rhoijres0=0;if (present(use_rhoijres)) use_rhoijres0=use_rhoijres
1130    if (paral_atom) then
1131      call pawrhoij_alloc(pawrhoij,cplex,nspden_rhoij,nspinor,nsppol,typat,&
1132 &     ngrhoij=ngrhoij0,nlmnmix=nlmnmix0,use_rhoij_=use_rhoij_0,use_rhoijres=use_rhoijres0,&
1133 &     pawtab=pawtab,comm_atom=my_comm_atom,mpi_atmtab=my_atmtab)
1134    else
1135      call pawrhoij_alloc(pawrhoij,cplex,nspden_rhoij,nspinor,nsppol,typat,pawtab=pawtab,&
1136 &     ngrhoij=ngrhoij0,nlmnmix=nlmnmix0,use_rhoij_=use_rhoij_0,use_rhoijres=use_rhoijres0)
1137    end if
1138  end if
1139 
1140  do iatom_rhoij=1,my_natom
1141    iatom=iatom_rhoij;if (paral_atom) iatom=my_atmtab(iatom_rhoij)
1142    itypat=typat(iatom)
1143    nselect=0
1144    ABI_ALLOCATE(lnspinat,(pawtab(itypat)%basis_size))
1145    lnspinat=-1
1146 ! Determine occupancies of each orbital
1147    if (nspden_rhoij==2) then
1148      ABI_ALLOCATE(occ,(pawtab(itypat)%basis_size))
1149      occ=zero
1150      do jlmn=1,pawtab(itypat)%lmn_size
1151        ln=pawtab(itypat)%indlmn(5,jlmn)
1152        klmn=jlmn*(jlmn+1)/2
1153        occ(ln)=occ(ln)+pawtab(itypat)%rhoij0(klmn)
1154      end do
1155      do ln=1,pawtab(itypat)%basis_size
1156        if(pawtab(itypat)%orbitals(ln)==0.and.occ(ln)==1) lnspinat(ln)=ln
1157        if(pawtab(itypat)%orbitals(ln)==1.and.(occ(ln)>=1.and.occ(ln)<=5)) lnspinat(ln)=ln
1158        if(pawtab(itypat)%orbitals(ln)==2.and.(occ(ln)>=1.and.occ(ln)<=9)) lnspinat(ln)=ln
1159        if(pawtab(itypat)%orbitals(ln)==3.and.(occ(ln)>=1.and.occ(ln)<=13)) lnspinat(ln)=ln
1160      end do
1161      ABI_DEALLOCATE(occ)
1162    end if
1163    lnspinat0=maxval(lnspinat)
1164    lnspinat0=-1
1165 
1166 !  Determine Z (trace of rhoij0 or part of it)
1167    zz=zero
1168    do jlmn=1,pawtab(itypat)%lmn_size
1169      jl=pawtab(itypat)%indlmn(1,jlmn)
1170      ln=pawtab(itypat)%indlmn(5,jlmn)
1171      j0lmn=jlmn*(jlmn-1)/2
1172      test_lnspinat=(lnspinat0==-1.or.lnspinat(ln)==ln)
1173      test_pawu=(lpawu(itypat)==-1.or.lpawu(itypat)==jl)
1174      test_exexch=(lexexch(itypat)==-1.or.lexexch(itypat)==jl)
1175      do ilmn=1,jlmn
1176        klmn=j0lmn+ilmn
1177        if ((ilmn==jlmn).and.test_pawu.and.test_exexch.and.test_lnspinat) &
1178 &       zz=zz+pawtab(itypat)%rhoij0(klmn)
1179      end do
1180    end do
1181 
1182 !  Compute rhoij from tabulated value and magnetization
1183    do ispden=1,nspden_rhoij
1184 
1185      zratio=zero
1186      roshift=one
1187      ratio=one
1188      if (nspden_rhoij==2) then
1189        ratio=half
1190        if ((spinat(3,iatom)>zero.and.ispden==1).or.&
1191 &       (spinat(3,iatom)<zero.and.ispden==2)) then
1192          if(abs(zz)>tol12)then
1193            zratio=two*abs(spinat(3,iatom))/zz
1194          else
1195            zratio=zero
1196          end if
1197        end if
1198      else if (nspden_rhoij==4.and.ispden>=2) then
1199        roshift=zero
1200        if(abs(zz)>tol12)then
1201          zratio=spinat(ispden-1,iatom)/zz
1202        else
1203          zratio=zero
1204        end if
1205      end if
1206 
1207      nselect=0;nselect1=1-cplex
1208      do jlmn=1,pawtab(itypat)%lmn_size
1209        jl=pawtab(itypat)%indlmn(1,jlmn)
1210        ln=pawtab(itypat)%indlmn(5,jlmn)
1211        j0lmn=jlmn*(jlmn-1)/2
1212        test_lnspinat=(lnspinat0==-1.or.lnspinat(ln)==ln)
1213        test_pawu=(lpawu(itypat)==-1.or.lpawu(itypat)==jl)
1214        test_exexch=(lexexch(itypat)==-1.or.lexexch(itypat)==jl)
1215        do ilmn=1,jlmn
1216          klmn=j0lmn+ilmn
1217          ro=pawtab(itypat)%rhoij0(klmn)
1218          if ((ilmn==jlmn).and.test_pawu.and.test_exexch.and.test_lnspinat) then
1219            ro=ro*ratio*(roshift+zratio)
1220          else
1221            ro=ro*ratio*roshift
1222          end if
1223 
1224          klmn1=cplex*(klmn-1)+1
1225          if (abs(ro)>tol10) then
1226            pawrhoij(iatom_rhoij)%rhoijp(klmn1,ispden)=ro
1227          else
1228            pawrhoij(iatom_rhoij)%rhoijp(klmn1,ispden)=zero
1229          end if
1230 
1231          if (ispden==nspden_rhoij) then
1232            if (any(abs(pawrhoij(iatom_rhoij)%rhoijp(klmn1,:))>tol10)) then
1233              nselect=nselect+1;nselect1=nselect1+cplex
1234              pawrhoij(iatom_rhoij)%rhoijselect(nselect)=klmn
1235              do jspden=1,nspden_rhoij
1236                pawrhoij(iatom_rhoij)%rhoijp(nselect1,jspden)=pawrhoij(iatom_rhoij)%rhoijp(klmn1,jspden)
1237              end do
1238            end if
1239          end if
1240 
1241        end do
1242      end do
1243 
1244    end do
1245    pawrhoij(iatom_rhoij)%nrhoijsel=nselect
1246 
1247 !  Non-collinear magnetism: avoid zero magnetization, because it produces numerical instabilities
1248 !    Add a small real to the magnetization ; not yet activated => must be tested.
1249 !   if (pawrhoij(iatom_rhoij)%nspden==4.and.spinat_zero) then
1250 !     pawrhoij(iatom_rhoij)%rhoijp(:,4)=pawrhoij(iatom_rhoij)%rhoijp(:,4)+tol10
1251 !   end if
1252    ABI_DEALLOCATE(lnspinat)
1253  end do ! iatom_rhoij
1254 
1255 !Destroy atom table used for parallelism
1256  call free_my_atmtab(my_atmtab,my_atmtab_allocated)
1257 
1258  DBG_EXIT("COLL")
1259 
1260 end subroutine initrhoij

m_paw_occupancies/m_paw_occupancies [ Modules ]

[ Top ] [ Modules ]

NAME

  m_paw_occupancies

FUNCTION

  This module contains routines related to the computation of PAW on-site occupancies (rhoij).

COPYRIGHT

 Copyright (C) 2018-2018 ABINIT group (FJ, MT)
 This file is distributed under the terms of the
 GNU General Public License, see ~abinit/COPYING
 or http://www.gnu.org/copyleft/gpl.txt .

SOURCE

17 #if defined HAVE_CONFIG_H
18 #include "config.h"
19 #endif
20 
21 #include "abi_common.h"
22 
23 MODULE m_paw_occupancies
24 
25  use defs_basis
26  use defs_abitypes
27  use m_abicore
28  use m_errors
29  use m_xmpi
30 
31  use m_pawtab,     only : pawtab_type
32  use m_pawrhoij,   only : pawrhoij_type,pawrhoij_init_unpacked,pawrhoij_mpisum_unpacked, &
33 &                         pawrhoij_alloc,pawrhoij_free,pawrhoij_get_nspden
34  use m_pawcprj,    only : pawcprj_type,pawcprj_alloc,pawcprj_get, &
35 &                         pawcprj_gather_spin, pawcprj_free
36  use m_paw_io,     only : pawio_print_ij
37  use m_paral_atom, only : get_my_atmtab,free_my_atmtab
38  use m_paw_dmft,   only : paw_dmft_type
39  use m_mpinfo,     only : proc_distrb_cycle
40 
41  implicit none
42 
43  private
44 
45 !public procedures.
46  public :: pawmkrhoij  ! Compute the PAW occupancies rhoij
47  public :: pawaccrhoij ! Accumulate the contribution of one band to the PAW occupancies rhoij
48  public :: initrhoij   ! Initialize the PAW occupancies rhoij from atomic data
49 
50 CONTAINS  !========================================================================================

m_paw_occupancies/pawaccrhoij [ Functions ]

[ Top ] [ m_paw_occupancies ] [ Functions ]

NAME

 pawaccrhoij

FUNCTION

 Accumulate the PAW quantities rhoij (augmentation occupancies)
 or their 1-st order change or their gradient vs r
 Add the contribution of a given k-point and band
 Remember: for each atom, rho_ij=Sum_{n,k} {occ(n,k)*<Cnk|p_i><p_j|Cnk>}

INPUTS

  atindx(natom)=index table for atoms (sorted-->random), inverse of atindx.
  cplex: if 1, WFs (or 1st-order WFs) are REAL, if 2, COMPLEX
  cwaveprj(natom,nspinor)= LEFT wave function at given n,k
                         projected with non-local projectors: cwaveprj=<p_i|Cnk>
  cwaveprj1(natom,nspinor)= RIGHT wave function at n,k,q
                          projected with non-local projectors: cwaveprj1=<p_i|C1nk,q>
                          * USED for RF  : C1nk is the first-order wave function
                          * USED for DMFT: C1nk is the RIGHT wave function
                          * NOT USED in usual GS case; can be set to cwaveprj in that case
  ipert=index of perturbation (RF only, i.e. option=2)
  isppol=index of current spin component
  mpi_atmtab(:)=--optional-- indexes of the atoms treated by current proc
  comm_atom=--optional-- MPI communicator over atoms
  my_natom=number of atoms treated by current processor
  natom=number of atoms in cell
  nspinor=number of spinorial components (on current proc)
  occ_k=occupation number for current band n,k
  option: choice of calculation:
          1: update rhoij (Ground-State)
          2: update 1-st order rhoij (Response Function) according to ipert
          3: update gradients of rhoij with respect to r,strain of both
  usetimerev=.TRUE. if time-reversal symmetry is used (WF(-k)=Conjg[WF(k)])
  wtk_k=weight assigned to current k-point

SIDE EFFECTS

  pawrhoij(natom) <type(pawrhoij_type)>= GS: paw rhoij occupancies and related data
                                         RF: 1-st order paw rhoij occupancies and related data
  On output, has been updated with the contribution of current n,k
    === option=1:
        pawrhoij(:)%rhoij_(lmn2_size,nspden)=      (non symetrized)
            Sum_{n,k} {occ(n,k)*conjugate[cprj_nk(ii)].cprj_nk(jj)}
    === option=2:
        pawrhoij(:)%rhoij_(lmn2_size,nspden)=      (non symetrized)
            Sum_{n,k} {occ(n,k)*(conjugate[cprj_nk(ii)].cprj1_nk,q(jj)
                                 conjugate[cprj_nk(jj)].cprj1_nk,q(ii)}
          + Sum_{n,k} {occ(n,k)*(conjugate[dcprj_nk(ii)/dlambda].cprj_nk(jj)
                                +conjugate[cprj_nk(ii)].dcprj_nk(jj)/dlambda)}
    === option=3:
        pawrhoij(:)%grhoij(lmn2_size,mu,nspden)=   (non symetrized)
            Sum_{n,k} {occ(n,k)*(conjugate[dcprj_nk(ii)/dr_mu].cprj_nk(jj)
                                +conjugate[cprj_nk(ii)].dcprj_nk(jj)/dr_mu)}

PARENTS

      d2frnl,dfpt_accrho,energy,pawmkrhoij,posdoppler,wfd_pawrhoij

CHILDREN

      free_my_atmtab,get_my_atmtab

SOURCE

 493  subroutine pawaccrhoij(atindx,cplex,cwaveprj,cwaveprj1,ipert,isppol,my_natom,natom,&
 494 &                       nspinor,occ_k,option,pawrhoij,usetimerev,wtk_k,occ_k_2, &
 495 &                       comm_atom,mpi_atmtab ) ! optional (parallelism)
 496 
 497 
 498 !This section has been created automatically by the script Abilint (TD).
 499 !Do not modify the following lines by hand.
 500 #undef ABI_FUNC
 501 #define ABI_FUNC 'pawaccrhoij'
 502 !End of the abilint section
 503 
 504  implicit none
 505 
 506 !Arguments ---------------------------------------------
 507 !scalars
 508  integer,intent(in) :: cplex,ipert,isppol,my_natom,natom,nspinor,option
 509  integer,optional,intent(in) :: comm_atom
 510  logical,intent(in) :: usetimerev
 511  real(dp),intent(in) :: occ_k,wtk_k
 512  real(dp),optional,intent(in) :: occ_k_2
 513 !arrays
 514  integer,intent(in) :: atindx(natom)
 515  integer,optional,target,intent(in) :: mpi_atmtab(:)
 516  type(pawcprj_type),intent(in) :: cwaveprj(natom,nspinor),cwaveprj1(natom,nspinor)
 517  type(pawrhoij_type),intent(inout) :: pawrhoij(my_natom)
 518 
 519 !Local variables ---------------------------------------
 520 !scalars
 521  integer :: cplex_rhoij,iatm,iatom,iatom1,ilmn,iplex,j0lmn,jlmn,klmn,klmn_im,klmn_re
 522  integer :: mu,my_comm_atom,ncpgr,nspden_rhoij
 523  logical :: compute_impart,compute_impart_cplex,substract_diagonal
 524  logical :: my_atmtab_allocated,paral_atom
 525  real(dp) :: ro11_im,ro11_re,ro12_im,ro12_re,ro21_im,ro21_re,ro22_im,ro22_re,weight,weight_2
 526  character(len=500) :: message
 527 !arrays
 528  integer,pointer :: my_atmtab(:)
 529  real(dp) :: cpi0(2,nspinor),cpi1(2,nspinor),cpj0(2,nspinor),cpj1(2,nspinor)
 530  real(dp) :: dcpi0(2,nspinor,9),dcpj0(2,nspinor,9)
 531 
 532 ! ***********************************************************************
 533 
 534  DBG_ENTER("COLL")
 535 
 536  if (my_natom==0) return
 537 
 538  ncpgr=0
 539  if (option==2.and.(ipert<=natom.or.ipert==natom+3.or.ipert==natom+4)) ncpgr=1
 540  if (option==3) ncpgr=cwaveprj(1,1)%ncpgr
 541 !Tests
 542  if(option==2.and.(ipert==natom+1.or.ipert==natom+10.or.ipert==natom+11)) then
 543    message = ' not relevant for ipert=natom+1 or ipert=natom+10 or ipert=natom+11 !'
 544    MSG_BUG(message)
 545  end if
 546  if(option==2.and.cwaveprj(1,1)%ncpgr<ncpgr) then
 547    message = ' Error on cwaveprj1 factors derivatives !'
 548    MSG_BUG(message)
 549  end if
 550  if(option==3.and.cwaveprj(1,1)%ncpgr/=ncpgr) then
 551    message = ' Error on cwaveprj factors derivatives !'
 552    MSG_BUG(message)
 553  end if
 554 
 555 !Set up parallelism over atoms
 556  paral_atom=(present(comm_atom).and.(my_natom/=natom))
 557  nullify(my_atmtab);if (present(mpi_atmtab)) my_atmtab => mpi_atmtab
 558  my_comm_atom=xmpi_comm_self;if (present(comm_atom)) my_comm_atom=comm_atom
 559  call get_my_atmtab(my_comm_atom,my_atmtab,my_atmtab_allocated,paral_atom,natom,&
 560 & my_natom_ref=my_natom)
 561 
 562  weight=wtk_k*occ_k
 563  weight_2=zero;if(present(occ_k_2).and.nspinor==2) weight_2=wtk_k*occ_k_2
 564  if (pawrhoij(1)%nspden==2.and.pawrhoij(1)%nsppol==1.and.nspinor==1) weight=half*weight
 565  if (pawrhoij(1)%nspden==2.and.pawrhoij(1)%nsppol==1.and.nspinor==1.and.present(occ_k_2)) weight_2=half*weight_2
 566 
 567  if (option==1) then
 568 
 569 !  ==================================================================
 570 !  === OPTION 1: Accumulate (n,k) contribution to rhoij =============
 571 !  ==================================================================
 572    compute_impart=((.not.usetimerev).and.(pawrhoij(1)%cplex==2))
 573    compute_impart_cplex=((compute_impart).and.(cplex==2))
 574    if (nspinor==1) then
 575      cplex_rhoij=pawrhoij(1)%cplex
 576      if (cplex_rhoij==1) then
 577        do iatom=1,my_natom
 578          iatom1=iatom;if (paral_atom) iatom1=my_atmtab(iatom)
 579          iatm=atindx(iatom1)
 580          do jlmn=1,pawrhoij(iatom)%lmn_size
 581            j0lmn=jlmn*(jlmn-1)/2
 582            cpj0(1:cplex,1)=cwaveprj(iatm,1)%cp(1:cplex,jlmn)
 583            do ilmn=1,jlmn
 584              klmn=j0lmn+ilmn
 585              cpi0(1:cplex,1)=cwaveprj1(iatm,1)%cp(1:cplex,ilmn)
 586              ro11_re=zero
 587              do iplex=1,cplex
 588                ro11_re=ro11_re+cpi0(iplex,1)*cpj0(iplex,1)
 589              end do
 590              pawrhoij(iatom)%rhoij_(klmn,isppol)=pawrhoij(iatom)%rhoij_(klmn,isppol)+weight*ro11_re
 591            end do
 592          end do
 593        end do
 594      else
 595        do iatom=1,my_natom
 596          iatom1=iatom;if (paral_atom) iatom1=my_atmtab(iatom)
 597          iatm=atindx(iatom1)
 598          do jlmn=1,pawrhoij(iatom)%lmn_size
 599            j0lmn=jlmn*(jlmn-1)/2
 600            cpj0(1:cplex,1)=cwaveprj(iatm,1)%cp(1:cplex,jlmn)
 601            do ilmn=1,jlmn
 602              klmn=j0lmn+ilmn
 603              klmn_re=cplex_rhoij*(klmn-1)+1
 604              cpi0(1:cplex,1)=cwaveprj1(iatm,1)%cp(1:cplex,ilmn)
 605              ro11_re=zero
 606              do iplex=1,cplex
 607                ro11_re=ro11_re+cpi0(iplex,1)*cpj0(iplex,1)
 608              end do
 609              pawrhoij(iatom)%rhoij_(klmn_re,isppol)=pawrhoij(iatom)%rhoij_(klmn_re,isppol)+weight*ro11_re
 610              if (compute_impart_cplex) then
 611                klmn_im=klmn_re+1
 612                ro11_im=cpi0(1,1)*cpj0(2,1)-cpi0(2,1)*cpj0(1,1)
 613                pawrhoij(iatom)%rhoij_(klmn_im,isppol)=pawrhoij(iatom)%rhoij_(klmn_im,isppol)+weight*ro11_im
 614              end if
 615            end do
 616          end do
 617        end do
 618      end if
 619    else ! nspinor=2
 620      do iatom=1,my_natom
 621        iatom1=iatom;if (paral_atom) iatom1=my_atmtab(iatom)
 622        iatm=atindx(iatom1)
 623        cplex_rhoij=pawrhoij(iatom)%cplex
 624        nspden_rhoij=pawrhoij(iatom)%nspden
 625        do jlmn=1,pawrhoij(iatom)%lmn_size
 626          j0lmn=jlmn*(jlmn-1)/2
 627          cpj0(1:cplex,1)=cwaveprj(iatm,1)%cp(1:cplex,jlmn)
 628          cpj0(1:cplex,2)=cwaveprj(iatm,2)%cp(1:cplex,jlmn)
 629          do ilmn=1,jlmn
 630            klmn=j0lmn+ilmn
 631            klmn_re=cplex_rhoij*(klmn-1)+1
 632            cpi0(1:cplex,1)=cwaveprj1(iatm,1)%cp(1:cplex,ilmn)
 633            cpi0(1:cplex,2)=cwaveprj1(iatm,2)%cp(1:cplex,ilmn)
 634            ro11_re=zero;ro22_re=zero
 635            ro12_re=zero;ro21_re=zero
 636            ro12_im=zero;ro21_im=zero
 637            do iplex=1,cplex
 638              ro11_re=ro11_re+cpi0(iplex,1)*cpj0(iplex,1)
 639              ro22_re=ro22_re+cpi0(iplex,2)*cpj0(iplex,2)
 640            end do
 641            pawrhoij(iatom)%rhoij_(klmn_re,1)=pawrhoij(iatom)%rhoij_(klmn_re,1)+weight*(ro11_re+ro22_re)
 642            if (nspden_rhoij>1) then
 643              do iplex=1,cplex
 644                ro12_re=ro12_re+cpi0(iplex,2)*cpj0(iplex,1)
 645                ro21_re=ro21_re+cpi0(iplex,1)*cpj0(iplex,2)
 646              end do
 647              pawrhoij(iatom)%rhoij_(klmn_re,4)=pawrhoij(iatom)%rhoij_(klmn_re,4)+weight*(ro11_re-ro22_re)
 648              pawrhoij(iatom)%rhoij_(klmn_re,2)=pawrhoij(iatom)%rhoij_(klmn_re,2)+weight*(ro12_re+ro21_re)
 649              if (cplex==2) then
 650                ro12_im=cpi0(1,2)*cpj0(2,1)-cpi0(2,2)*cpj0(1,1)
 651                ro21_im=cpi0(1,1)*cpj0(2,2)-cpi0(2,1)*cpj0(1,2)
 652                pawrhoij(iatom)%rhoij_(klmn_re,3)=pawrhoij(iatom)%rhoij_(klmn_re,3)+weight*(ro21_im-ro12_im)
 653              end if
 654            end if
 655            if (compute_impart) then
 656              klmn_im=klmn_re+1
 657              if (nspden_rhoij>1) pawrhoij(iatom)%rhoij_(klmn_im,3)=pawrhoij(iatom)%rhoij_(klmn_im,3)+weight*(ro12_re-ro21_re)
 658              if (cplex==2) then
 659                ro11_im=cpi0(1,1)*cpj0(2,1)-cpi0(2,1)*cpj0(1,1)
 660                ro22_im=cpi0(1,2)*cpj0(2,2)-cpi0(2,2)*cpj0(1,2)
 661                pawrhoij(iatom)%rhoij_(klmn_im,1)=pawrhoij(iatom)%rhoij_(klmn_im,1)+weight*(ro11_im+ro22_im)
 662                if (nspden_rhoij>1) then
 663                  pawrhoij(iatom)%rhoij_(klmn_im,4)=pawrhoij(iatom)%rhoij_(klmn_im,4)+weight*(ro11_im-ro22_im)
 664                  pawrhoij(iatom)%rhoij_(klmn_im,2)=pawrhoij(iatom)%rhoij_(klmn_im,2)+weight*(ro12_im+ro21_im)
 665                end if
 666                if (present(occ_k_2).and.nspinor==2) then
 667                  pawrhoij(iatom)%rhoij_(klmn_re,1)=pawrhoij(iatom)%rhoij_(klmn_re,1)+weight_2*(-ro11_im-ro22_im)
 668                  pawrhoij(iatom)%rhoij_(klmn_im,1)=pawrhoij(iatom)%rhoij_(klmn_im,1)+weight_2*( ro11_re+ro22_re)
 669                  pawrhoij(iatom)%rhoij_(klmn_re,2)=pawrhoij(iatom)%rhoij_(klmn_re,2)+weight_2*(-ro21_im-ro12_im)
 670                  pawrhoij(iatom)%rhoij_(klmn_im,2)=pawrhoij(iatom)%rhoij_(klmn_im,2)+weight_2*( ro21_re+ro12_re)
 671                  pawrhoij(iatom)%rhoij_(klmn_re,3)=pawrhoij(iatom)%rhoij_(klmn_re,3)+weight_2*(-ro12_re+ro21_re)
 672                  pawrhoij(iatom)%rhoij_(klmn_im,3)=pawrhoij(iatom)%rhoij_(klmn_im,3)+weight_2*(-ro12_im+ro21_im)
 673                  pawrhoij(iatom)%rhoij_(klmn_re,4)=pawrhoij(iatom)%rhoij_(klmn_re,4)+weight_2*(-ro11_im+ro22_im)
 674                  pawrhoij(iatom)%rhoij_(klmn_im,4)=pawrhoij(iatom)%rhoij_(klmn_im,4)+weight_2*( ro11_re-ro22_re)
 675                end if
 676              end if
 677            end if
 678          end do
 679        end do
 680      end do
 681    end if
 682 
 683  else if (option==2) then
 684 
 685 !  ==================================================================
 686 !  === OPTION 2: Accumulate (n,k) contribution to 1st-order rhoij ===
 687 !  ==================================================================
 688 
 689    compute_impart=(pawrhoij(1)%cplex==2)
 690    compute_impart_cplex=((pawrhoij(1)%cplex==2).and.(cplex==2))
 691    substract_diagonal=(ipert==natom+3)
 692 
 693    if (compute_impart_cplex) then
 694      if (.not.allocated(pawrhoij(1)%rhoijim)) then
 695        MSG_BUG("pawrhoij(:)%rhoijim must be allocated!")
 696      end if
 697    end if
 698 
 699 !  Accumulate (n,k) contribution to rhoij1
 700 !  due to derivative of wave-function
 701    if (nspinor==1) then
 702      do iatom=1,my_natom
 703        iatom1=iatom;if (paral_atom) iatom1=my_atmtab(iatom)
 704        iatm=atindx(iatom1)
 705        cplex_rhoij=pawrhoij(iatom)%cplex
 706        do jlmn=1,pawrhoij(iatom)%lmn_size
 707          j0lmn=jlmn*(jlmn-1)/2
 708          cpj0(1:2,1)=cwaveprj (iatm,1)%cp(1:2,jlmn)
 709          cpj1(1:2,1)=cwaveprj1(iatm,1)%cp(1:2,jlmn)
 710          do ilmn=1,jlmn
 711            klmn=j0lmn+ilmn
 712            klmn_re=cplex_rhoij*(klmn-1)+1
 713            cpi0(1:2,1)=cwaveprj (iatm,1)%cp(1:2,ilmn)
 714            cpi1(1:2,1)=cwaveprj1(iatm,1)%cp(1:2,ilmn)
 715            ro11_re=zero
 716            do iplex=1,cplex
 717              ro11_re=ro11_re+cpi0(iplex,1)*cpj1(iplex,1)+cpj0(iplex,1)*cpi1(iplex,1)
 718            end do
 719            pawrhoij(iatom)%rhoij_(klmn_re,isppol)=pawrhoij(iatom)%rhoij_(klmn_re,isppol)+weight*ro11_re
 720            if (compute_impart_cplex) then
 721              klmn_im=klmn_re+1
 722              ro11_im=cpi0(1,1)*cpj1(2,1)-cpi0(2,1)*cpj1(1,1)+cpj0(1,1)*cpi1(2,1)-cpj0(2,1)*cpi1(1,1)
 723              pawrhoij(iatom)%rhoij_(klmn_im,isppol)=pawrhoij(iatom)%rhoij_(klmn_im,isppol)+weight*ro11_im
 724            end if
 725          end do
 726        end do
 727      end do
 728    else ! nspinor=2
 729      do iatom=1,my_natom
 730        iatom1=iatom;if (paral_atom) iatom1=my_atmtab(iatom)
 731        iatm=atindx(iatom1)
 732        cplex_rhoij=pawrhoij(iatom)%cplex
 733        nspden_rhoij=pawrhoij(iatom)%nspden
 734        do jlmn=1,pawrhoij(iatom)%lmn_size
 735          j0lmn=jlmn*(jlmn-1)/2
 736          cpj0(1:2,1)=cwaveprj (iatm,1)%cp(1:2,jlmn)
 737          cpj0(1:2,2)=cwaveprj (iatm,2)%cp(1:2,jlmn)
 738          cpj1(1:2,1)=cwaveprj1(iatm,1)%cp(1:2,jlmn)
 739          cpj1(1:2,2)=cwaveprj1(iatm,2)%cp(1:2,jlmn)
 740          do ilmn=1,jlmn
 741            klmn=j0lmn+ilmn
 742            klmn_re=cplex_rhoij*(klmn-1)+1
 743            cpi0(1:2,1)=cwaveprj (iatm,1)%cp(1:2,ilmn)
 744            cpi0(1:2,2)=cwaveprj (iatm,2)%cp(1:2,ilmn)
 745            cpi1(1:2,1)=cwaveprj1(iatm,1)%cp(1:2,ilmn)
 746            cpi1(1:2,2)=cwaveprj1(iatm,2)%cp(1:2,ilmn)
 747            ro11_re=zero;ro22_re=zero
 748            ro12_re=zero;ro21_re=zero
 749            ro12_im=zero;ro21_im=zero
 750            do iplex=1,cplex
 751              ro11_re=ro11_re+cpj0(iplex,1)*cpi1(iplex,1)+cpi0(iplex,1)*cpj1(iplex,1)
 752              ro22_re=ro22_re+cpj0(iplex,2)*cpi1(iplex,2)+cpi0(iplex,2)*cpj1(iplex,2)
 753            end do
 754            pawrhoij(iatom)%rhoij_(klmn_re,1)=pawrhoij(iatom)%rhoij_(klmn_re,1)+weight*(ro11_re+ro22_re)
 755            if (nspden_rhoij>1) then
 756              do iplex=1,cplex
 757                ro12_re=ro12_re+cpj0(iplex,1)*cpi1(iplex,2)+cpi0(iplex,2)*cpj1(iplex,1)
 758                ro21_re=ro21_re+cpj0(iplex,2)*cpi1(iplex,1)+cpi0(iplex,1)*cpj1(iplex,2)
 759              end do
 760              pawrhoij(iatom)%rhoij_(klmn_re,4)=pawrhoij(iatom)%rhoij_(klmn_re,4)+weight*(ro11_re-ro22_re)
 761              pawrhoij(iatom)%rhoij_(klmn_re,2)=pawrhoij(iatom)%rhoij_(klmn_re,2)+weight*(ro12_re+ro21_re)
 762              if (cplex==2) then
 763                ro12_im=cpj0(1,1)*cpi1(2,2)-cpi1(1,1)*cpj0(2,2)+cpi0(1,2)*cpj1(2,1)-cpj1(1,2)*cpi0(2,1)
 764                ro21_im=cpj0(1,2)*cpi1(2,1)-cpi1(1,2)*cpj0(2,1)+cpi0(1,1)*cpj1(2,2)-cpj1(1,1)*cpi0(2,2)
 765                pawrhoij(iatom)%rhoij_(klmn_re,3)=pawrhoij(iatom)%rhoij_(klmn_re,3)+weight*(ro21_im-ro12_im)
 766              end if
 767            end if
 768            if (compute_impart) then
 769              klmn_im=klmn_re+1
 770              if (nspden_rhoij>1) pawrhoij(iatom)%rhoij_(klmn_re,3)=pawrhoij(iatom)%rhoij_(klmn_re,3)+weight*(ro12_re-ro21_re)
 771              if (cplex==2) then
 772                ro11_im=cpj0(1,1)*cpi1(2,1)-cpi1(1,1)*cpj0(2,1)+cpi0(1,1)*cpj1(2,1)-cpj1(1,1)*cpi0(2,1)
 773                ro22_im=cpj0(1,2)*cpi1(2,2)-cpi1(1,2)*cpj0(2,2)+cpi0(1,2)*cpj1(2,2)-cpj1(1,2)*cpi0(2,2)
 774                pawrhoij(iatom)%rhoij_(klmn_im,1)=pawrhoij(iatom)%rhoij_(klmn_im,1)+weight*(ro11_im+ro22_im)
 775                if (nspden_rhoij>1) then
 776                  pawrhoij(iatom)%rhoij_(klmn_im,4)=pawrhoij(iatom)%rhoij_(klmn_im,4)+weight*(ro11_im-ro22_im)
 777                  pawrhoij(iatom)%rhoij_(klmn_re,2)=pawrhoij(iatom)%rhoij_(klmn_re,2)+weight*(ro12_im+ro21_im)
 778                end if
 779              end if
 780            end if
 781          end do
 782        end do
 783      end do
 784    end if
 785 
 786 !  Accumulate (n,k) contribution to rhoij1
 787 !  due to derivative of projectors
 788    if (ipert/=natom+2) then
 789      if (nspinor==1) then
 790        do iatom=1,my_natom
 791          iatom1=iatom;if (paral_atom) iatom1=my_atmtab(iatom)
 792          iatm=atindx(iatom1)
 793          if (ipert<=natom.and.iatom/=ipert) cycle
 794          cplex_rhoij=pawrhoij(iatom)%cplex
 795          do jlmn=1,pawrhoij(iatom)%lmn_size
 796            j0lmn=jlmn*(jlmn-1)/2
 797            cpj0 (1:2,1)  =cwaveprj(iatm,1)%cp (1:2  ,jlmn)
 798            dcpj0(1:2,1,1)=cwaveprj(iatm,1)%dcp(1:2,1,jlmn)
 799            do ilmn=1,jlmn
 800              klmn=j0lmn+ilmn
 801              klmn_re=cplex_rhoij*(klmn-1)+1
 802              cpi0 (1:2,1)  =cwaveprj(iatm,1)%cp (1:2  ,ilmn)
 803              dcpi0(1:2,1,1)=cwaveprj(iatm,1)%dcp(1:2,1,ilmn)
 804              ro11_re=zero
 805              do iplex=1,cplex
 806                ro11_re=ro11_re+dcpi0(iplex,1,1)*cpj0(iplex,1)+cpi0(iplex,1)*dcpj0(iplex,1,1)
 807              end do
 808              if (substract_diagonal) then
 809                do iplex=1,cplex
 810                  ro11_re=ro11_re-cpi0(iplex,1)*cpj0(iplex,1)
 811                end do
 812              end if
 813              pawrhoij(iatom)%rhoij_(klmn_re,isppol)=pawrhoij(iatom)%rhoij_(klmn_re,isppol)+weight*ro11_re
 814 !            This imaginary part does not have to be computed
 815 !            It is cancelled because rho_ij+rho_ji is stored in rho_ij
 816              if (compute_impart_cplex) then
 817                klmn_im=klmn_re+1
 818                ro11_im=dcpi0(1,1,1)*cpj0(2,1)-dcpi0(2,1,1)*cpj0(1,1)+cpi0(1,1)*dcpj0(2,1,1)-cpi0(2,1)*dcpj0(1,1,1)
 819                if (substract_diagonal) then
 820                  ro11_im=ro11_im-cpi0(1,1)*cpj0(2,1)+cpi0(2,1)*cpj0(1,1)
 821                end if
 822                pawrhoij(iatom)%rhoijim(klmn,isppol)=pawrhoij(iatom)%rhoijim(klmn,isppol)+weight*ro11_im
 823 !              pawrhoij(iatom)%rhoij_(klmn_im,isppol)=pawrhoij(iatom)%rhoij_(klmn_im,isppol)+weight*ro11_im
 824              end if
 825            end do
 826          end do
 827        end do
 828      else ! nspinor=2
 829        do iatom=1,my_natom
 830          iatom1=iatom;if (paral_atom) iatom1=my_atmtab(iatom)
 831          iatm=atindx(iatom1)
 832          if (ipert<=natom.and.iatom/=ipert) cycle
 833          cplex_rhoij=pawrhoij(iatom)%cplex
 834          nspden_rhoij=pawrhoij(iatom)%nspden
 835          do jlmn=1,pawrhoij(iatom)%lmn_size
 836            j0lmn=jlmn*(jlmn-1)/2
 837            cpj0 (1:2,1)  =cwaveprj(iatm,1)%cp (1:2  ,jlmn)
 838            dcpj0(1:2,1,1)=cwaveprj(iatm,1)%dcp(1:2,1,jlmn)
 839            cpj0 (1:2,2)  =cwaveprj(iatm,2)%cp (1:2  ,jlmn)
 840            dcpj0(1:2,2,1)=cwaveprj(iatm,2)%dcp(1:2,1,jlmn)
 841            do ilmn=1,jlmn
 842              klmn=j0lmn+ilmn
 843              klmn_re=cplex_rhoij*(klmn-1)+1
 844              cpi0 (1:2,1)  =cwaveprj(iatm,1)%cp (1:2  ,ilmn)
 845              dcpi0(1:2,1,1)=cwaveprj(iatm,1)%dcp(1:2,1,ilmn)
 846              cpi0 (1:2,2)  =cwaveprj(iatm,2)%cp (1:2  ,ilmn)
 847              dcpi0(1:2,2,1)=cwaveprj(iatm,2)%dcp(1:2,1,ilmn)
 848              ro11_re=zero;ro22_re=zero
 849              ro12_re=zero;ro21_re=zero
 850              ro12_im=zero;ro21_im=zero
 851              do iplex=1,cplex
 852                ro11_re=ro11_re+dcpi0(iplex,1,1)*cpj0(iplex,1)+cpi0(iplex,1)*dcpj0(iplex,1,1)
 853                ro22_re=ro22_re+dcpi0(iplex,2,1)*cpj0(iplex,2)+cpi0(iplex,2)*dcpj0(iplex,2,1)
 854              end do
 855              if (substract_diagonal) then
 856                do iplex=1,cplex
 857                  ro11_re=ro11_re-cpi0(iplex,1)*cpj0(iplex,1)
 858                  ro22_re=ro22_re-cpi0(iplex,2)*cpj0(iplex,2)
 859                end do
 860              end if
 861              pawrhoij(iatom)%rhoij_(klmn_re,1)=pawrhoij(iatom)%rhoij_(klmn_re,1)+weight*(ro11_re+ro22_re)
 862              if (nspden_rhoij>1) then
 863                do iplex=1,cplex
 864                  ro12_re=ro12_re+dcpi0(iplex,2,1)*cpj0(iplex,1)+cpi0(iplex,2)*dcpj0(iplex,1,1)
 865                  ro21_re=ro21_re+dcpi0(iplex,1,1)*cpj0(iplex,2)+cpi0(iplex,1)*dcpj0(iplex,2,1)
 866                end do
 867                if (substract_diagonal) then
 868                  do iplex=1,cplex
 869                    ro12_re=ro12_re-cpi0(iplex,2)*cpj0(iplex,1)
 870                    ro21_re=ro21_re-cpi0(iplex,1)*cpj0(iplex,2)
 871                  end do
 872                end if
 873                pawrhoij(iatom)%rhoij_(klmn_re,4)=pawrhoij(iatom)%rhoij_(klmn_re,4)+weight*(ro11_re-ro22_re)
 874                pawrhoij(iatom)%rhoij_(klmn_re,2)=pawrhoij(iatom)%rhoij_(klmn_re,2)+weight*(ro12_re+ro21_re)
 875                if (cplex==2) then
 876                  ro12_im=dcpi0(1,2,1)*cpj0(2,1)-dcpi0(2,2,1)*cpj0(1,1)+cpi0(1,2)*dcpj0(2,1,1)-cpi0(2,2)*dcpj0(1,1,1)
 877                  ro21_im=dcpi0(1,1,1)*cpj0(2,2)-dcpi0(2,1,1)*cpj0(1,2)+cpi0(1,1)*dcpj0(2,2,1)-cpi0(2,1)*dcpj0(1,2,1)
 878                  if (substract_diagonal) then
 879                    ro12_im=ro12_im-cpi0(1,2)*cpj0(2,1)+cpi0(2,2)*cpj0(1,1)
 880                    ro21_im=ro21_im-cpi0(1,1)*cpj0(2,2)+cpi0(2,1)*cpj0(1,2)
 881                  end if
 882                  pawrhoij(iatom)%rhoij_(klmn_re,3)=pawrhoij(iatom)%rhoij_(klmn_re,3)+weight*(ro21_im-ro12_im)
 883                end if
 884              end if
 885              if (compute_impart) then
 886                klmn_im=klmn_re+1
 887                if (nspden_rhoij>1) pawrhoij(iatom)%rhoij_(klmn_im,3)=pawrhoij(iatom)%rhoij_(klmn_im,3)+weight*(ro12_re-ro21_re)
 888                if (cplex==2) then
 889                  ro11_im=dcpi0(1,1,1)*cpj0(2,1)-dcpi0(2,1,1)*cpj0(1,1)+cpi0(1,1)*dcpj0(2,1,1)-cpi0(2,1)*dcpj0(1,1,1)
 890                  ro22_im=dcpi0(1,2,1)*cpj0(2,2)-dcpi0(2,2,1)*cpj0(1,2)+cpi0(1,2)*dcpj0(2,2,1)-cpi0(2,2)*dcpj0(1,2,1)
 891                  if (substract_diagonal) then
 892                    ro11_im=ro11_im-cpi0(1,1)*cpj0(2,1)+cpi0(2,1)*cpj0(1,1)
 893                    ro22_im=ro22_im-cpi0(1,2)*cpj0(2,2)+cpi0(2,2)*cpj0(1,2)
 894                  end if
 895                  pawrhoij(iatom)%rhoij_(klmn_im,1)=pawrhoij(iatom)%rhoij_(klmn_im,1)+weight*(ro11_im+ro22_im)
 896                  if (nspden_rhoij>1) then
 897                    pawrhoij(iatom)%rhoij_(klmn_im,4)=pawrhoij(iatom)%rhoij_(klmn_im,4)+weight*(ro11_im-ro22_im)
 898                    pawrhoij(iatom)%rhoij_(klmn_im,2)=pawrhoij(iatom)%rhoij_(klmn_im,2)+weight*(ro12_im+ro21_im)
 899                  end if
 900                end if
 901              end if
 902            end do
 903          end do
 904        end do
 905      end if
 906    end if
 907 
 908  else if (option==3) then
 909 
 910 !  ==================================================================
 911 !  === OPTION 3: Accumulate (n,k) contribution to drhoij/dr =========
 912 !  ==================================================================
 913 
 914    compute_impart=((.not.usetimerev).and.(pawrhoij(1)%cplex==2))
 915    compute_impart_cplex=((compute_impart).and.(cplex==2))
 916    if (nspinor==1) then
 917      do iatom=1,my_natom
 918        iatom1=iatom;if (paral_atom) iatom1=my_atmtab(iatom)
 919        iatm=atindx(iatom1)
 920        cplex_rhoij=pawrhoij(iatom)%cplex
 921        do jlmn=1,pawrhoij(iatom)%lmn_size
 922          j0lmn=jlmn*(jlmn-1)/2
 923          cpj0(1:cplex,1)         =cwaveprj(iatm,1)%cp (1:cplex,jlmn)
 924          dcpj0(1:cplex,1,1:ncpgr)=cwaveprj(iatm,1)%dcp(1:cplex,1:ncpgr,jlmn)
 925          do ilmn=1,jlmn
 926            klmn=j0lmn+ilmn
 927            klmn_re=cplex_rhoij*(klmn-1)+1
 928            cpi0(1:cplex,1)         =cwaveprj(iatm,1)%cp (1:cplex,ilmn)
 929            dcpi0(1:cplex,1,1:ncpgr)=cwaveprj(iatm,1)%dcp(1:cplex,1:ncpgr,ilmn)
 930            do mu=1,ncpgr
 931              ro11_re=zero
 932              do iplex=1,cplex
 933                ro11_re=ro11_re+dcpi0(iplex,1,mu)*cpj0(iplex,1)+cpi0(iplex,1)*dcpj0(iplex,1,mu)
 934              end do
 935              pawrhoij(iatom)%grhoij(mu,klmn_re,isppol)=pawrhoij(iatom)%grhoij(mu,klmn_re,isppol)+weight*ro11_re
 936            end do
 937            if (compute_impart_cplex) then
 938              klmn_im=klmn_re+1
 939              do mu=1,ncpgr
 940                ro11_im=dcpi0(1,1,mu)*cpj0(2,1)+cpi0(1,1)*dcpj0(2,1,mu)-dcpi0(2,1,mu)*cpj0(1,1)-cpi0(2,1)*dcpj0(1,1,mu)
 941                pawrhoij(iatom)%grhoij(mu,klmn_im,isppol)=pawrhoij(iatom)%grhoij(mu,klmn_im,isppol)+weight*ro11_im
 942              end do
 943            end if
 944          end do
 945        end do
 946      end do
 947    else ! nspinor=2
 948      do iatom=1,my_natom
 949        iatom1=iatom;if (paral_atom) iatom1=my_atmtab(iatom)
 950        iatm=atindx(iatom1)
 951        cplex_rhoij=pawrhoij(iatom)%cplex
 952        nspden_rhoij=pawrhoij(iatom)%nspden
 953        do jlmn=1,pawrhoij(iatom)%lmn_size
 954          j0lmn=jlmn*(jlmn-1)/2
 955          cpj0(1:cplex,1)     =cwaveprj(iatm,1)%cp (1:cplex,jlmn)
 956          cpj0(1:cplex,2)     =cwaveprj(iatm,2)%cp (1:cplex,jlmn)
 957          dcpj0(1:cplex,1,1:ncpgr)=cwaveprj(iatm,1)%dcp(1:cplex,1:ncpgr,jlmn)
 958          dcpj0(1:cplex,2,1:ncpgr)=cwaveprj(iatm,2)%dcp(1:cplex,1:ncpgr,jlmn)
 959          do ilmn=1,jlmn
 960            klmn=j0lmn+ilmn
 961            klmn_re=cplex_rhoij*(klmn-1)+1
 962            klmn_im=klmn_re+1
 963            cpi0(1:cplex,1)     =cwaveprj(iatm,1)%cp (1:cplex,ilmn)
 964            cpi0(1:cplex,2)     =cwaveprj(iatm,2)%cp (1:cplex,ilmn)
 965            dcpi0(1:cplex,1,1:ncpgr)=cwaveprj(iatm,1)%dcp(1:cplex,1:ncpgr,ilmn)
 966            dcpi0(1:cplex,2,1:ncpgr)=cwaveprj(iatm,2)%dcp(1:cplex,1:ncpgr,ilmn)
 967            do mu=1,ncpgr
 968              ro11_re=zero;ro22_re=zero
 969              ro12_re=zero;ro21_re=zero
 970              ro12_im=zero;ro21_im=zero
 971              do iplex=1,cplex
 972                ro11_re=ro11_re+dcpi0(iplex,1,mu)*cpj0(iplex,1)+cpi0(iplex,1)*dcpj0(iplex,1,mu)
 973                ro22_re=ro22_re+dcpi0(iplex,2,mu)*cpj0(iplex,2)+cpi0(iplex,2)*dcpj0(iplex,2,mu)
 974              end do
 975              pawrhoij(iatom)%grhoij(mu,klmn_re,1)=pawrhoij(iatom)%grhoij(mu,klmn_re,1)+weight*(ro11_re+ro22_re)
 976              if (nspden_rhoij>1) then
 977                do iplex=1,cplex
 978                  ro12_re=ro12_re+dcpi0(iplex,2,mu)*cpj0(iplex,1)+cpi0(iplex,2)*dcpj0(iplex,1,mu)
 979                  ro21_re=ro21_re+dcpi0(iplex,1,mu)*cpj0(iplex,2)+cpi0(iplex,1)*dcpj0(iplex,2,mu)
 980                end do
 981                pawrhoij(iatom)%grhoij(mu,klmn_re,4)=pawrhoij(iatom)%grhoij(mu,klmn_re,4)+weight*(ro11_re-ro22_re)
 982                pawrhoij(iatom)%grhoij(mu,klmn_re,2)=pawrhoij(iatom)%grhoij(mu,klmn_re,2)+weight*(ro12_re+ro21_re)
 983                if (cplex==2) then
 984                  ro12_im=dcpi0(1,2,mu)*cpj0(2,1)+cpi0(1,2)*dcpj0(2,1,mu)-dcpi0(2,2,mu)*cpj0(1,1)-cpi0(2,2)*dcpj0(1,1,mu)
 985                  ro21_im=dcpi0(1,1,mu)*cpj0(2,2)+cpi0(1,1)*dcpj0(2,2,mu)-dcpi0(2,1,mu)*cpj0(1,2)-cpi0(2,1)*dcpj0(1,2,mu)
 986                  pawrhoij(iatom)%grhoij(mu,klmn_re,3)=pawrhoij(iatom)%grhoij(mu,klmn_re,3)+weight*(ro21_im-ro12_im)
 987                end if
 988              end if
 989              if (compute_impart) then
 990                if (nspden_rhoij>1) then
 991                  pawrhoij(iatom)%grhoij(mu,klmn_im,3)=pawrhoij(iatom)%grhoij(mu,klmn_im,3)+weight*(ro12_re-ro21_re)
 992                end if
 993                if (cplex==2) then
 994                  ro11_im=dcpi0(1,1,mu)*cpj0(2,1)+cpi0(1,1)*dcpj0(2,1,mu)-dcpi0(2,1,mu)*cpj0(1,1)-cpi0(2,1)*dcpj0(1,1,mu)
 995                  ro22_im=dcpi0(1,2,mu)*cpj0(2,2)+cpi0(1,2)*dcpj0(2,2,mu)-dcpi0(2,2,mu)*cpj0(1,2)-cpi0(2,2)*dcpj0(1,2,mu)
 996                  pawrhoij(iatom)%grhoij(mu,klmn_im,1)=pawrhoij(iatom)%grhoij(mu,klmn_im,1)+weight*(ro11_im+ro22_im)
 997                  if (nspden_rhoij>1) then
 998                    pawrhoij(iatom)%grhoij(mu,klmn_im,4)=pawrhoij(iatom)%grhoij(mu,klmn_im,4)+weight*(ro11_im-ro22_im)
 999                    pawrhoij(iatom)%grhoij(mu,klmn_im,2)=pawrhoij(iatom)%grhoij(mu,klmn_im,2)+weight*(ro12_im+ro21_im)
1000                  end if
1001                end if
1002              end if
1003            end do
1004          end do
1005        end do
1006      end do
1007    end if
1008 
1009 !  End
1010  end if ! option
1011 
1012 !Destroy atom table used for parallelism
1013  call free_my_atmtab(my_atmtab,my_atmtab_allocated)
1014 
1015  DBG_EXIT("COLL")
1016 
1017 end subroutine pawaccrhoij

m_paw_occupancies/pawmkrhoij [ Functions ]

[ Top ] [ m_paw_occupancies ] [ Functions ]

NAME

 pawmkrhoij

FUNCTION

 Calculate the PAW quantities rhoij (augmentation occupancies)
 Remember:for each atom, rho_ij=Sum_{n,k} {occ(n,k)*<Cnk|p_i><p_j|Cnk>}

INPUTS

  atindx1(natom)=index table for atoms, inverse of atindx
  cprj(natom,mcprj)= wave functions projected with non-local projectors:
                     cprj_nk(i)=<p_i|Cnk> where p_i is a non-local projector.
  dimcprj(natom)=array of dimensions of array cprj (ordered by atom-type)
  istwfk(nkpt)=parameter that describes the storage of wfs
  kptopt=option for the generation of k points
  mband=maximum number of bands
  mband_cprj=maximum number of bands used in the dimensioning of cprj array (usually mband/nproc_band)
  mcprj=size of projected wave-functions array (cprj) =nspinor*mband*mkmem*nsppol
  mkmem =number of k points treated by this node.
  mpi_enreg=information about MPI parallelization
  natom=number of atoms in cell
  nband=number of bands for all k points
  nkpt=number of k points
  nspinor=number of spinorial components of the wavefunctions
  nsppol=1 for unpolarized, 2 for spin-polarized
  occ(mband*nkpt*nsppol)=occupation number for each band for each k
  paral_kgb=Flag related to the kpoint-band-fft parallelism
  paw_dmft  <type(paw_dmft_type)>= paw+dmft related data
  pawprtvol=control print volume and debugging output for PAW
  unpaw=unit number for cprj PAW data (if used)
  wtk(nkpt)=weight assigned to each k point

SIDE EFFECTS

  pawrhoij(natom) <type(pawrhoij_type)>= paw rhoij occupancies and related data
  On input: arrays dimensions
  On output:
    pawrhoij(:)%rhoij_(lmn2_size,nspden)=
          Sum_{n,k} {occ(n,k)*conjugate[cprj_nk(ii)].cprj_nk(jj)} (non symetrized)

PARENTS

      afterscfloop,scfcv,vtorho

CHILDREN

      pawaccrhoij,pawcprj_alloc,pawcprj_free,pawcprj_gather_spin,pawcprj_get
      pawio_print_ij,pawrhoij_free,pawrhoij_init_unpacked
      pawrhoij_mpisum_unpacked,wrtout

NOTES

  The cprj are distributed over band processors.
  Only the mod((iband-1)/mpi_enreg%bandpp,mpi_enreg%nproc_band) projectors
  are stored on each proc.

SOURCE

110  subroutine pawmkrhoij(atindx,atindx1,cprj,dimcprj,istwfk,kptopt,mband,mband_cprj,mcprj,mkmem,mpi_enreg,&
111 &                      natom,nband,nkpt,nspinor,nsppol,occ,paral_kgb,paw_dmft,&
112 &                      pawprtvol,pawrhoij,unpaw,usewvl,wtk)
113 
114 
115 !This section has been created automatically by the script Abilint (TD).
116 !Do not modify the following lines by hand.
117 #undef ABI_FUNC
118 #define ABI_FUNC 'pawmkrhoij'
119 !End of the abilint section
120 
121  implicit none
122 
123 !Arguments ---------------------------------------------
124 !scalars
125  integer,intent(in) :: kptopt,mband,mband_cprj,mcprj,mkmem,natom,nkpt,nspinor,nsppol
126  integer,intent(in) :: paral_kgb,pawprtvol,unpaw,usewvl
127  type(MPI_type),intent(in) :: mpi_enreg
128 !arrays
129  integer,intent(in) :: atindx(natom),atindx1(natom),dimcprj(natom),istwfk(nkpt)
130  integer,intent(in) :: nband(nkpt*nsppol)
131  real(dp),intent(in) :: occ(mband*nkpt*nsppol),wtk(nkpt)
132  type(pawcprj_type),target,intent(in) :: cprj(natom,mcprj)
133  type(paw_dmft_type),intent(in) :: paw_dmft
134  type(pawrhoij_type),intent(inout),target:: pawrhoij(:)
135 
136 !Local variables ---------------------------------------
137 !scalars
138  integer,parameter :: max_nband_cprj=100
139  integer :: bdtot_index,cplex
140  integer :: iatom,iatom_tot,ib,ib1,iband,iband1,ibc1,ibg,ib_this_proc,ierr
141  integer :: ikpt,iorder_cprj,isppol,jb_this_proc,jbg,me,my_nspinor,nband_k,nband_k_cprj
142  integer :: nbandc1,nband_k_cprj_read,nband_k_cprj_used,nprocband,nrhoij,nsp2
143  integer :: option,spaceComm,use_nondiag_occup_dmft
144  logical :: locc_test,paral_atom,usetimerev
145  real(dp) :: wtk_k
146  character(len=4) :: wrt_mode
147  character(len=500) :: msg
148 
149 !arrays
150  integer :: idum(0)
151  real(dp) :: occup(2)
152  character(len=8),parameter :: dspin(6)=(/"up      ","down    ","dens (n)","magn (x)","magn (y)","magn (z)"/)
153  type(pawcprj_type),allocatable :: cprj_tmp(:,:),cwaveprj(:,:),cwaveprjb(:,:)
154  type(pawcprj_type),pointer :: cprj_ptr(:,:)
155  type(pawrhoij_type),pointer :: pawrhoij_all(:)
156 
157 !************************************************************************
158 
159  DBG_ENTER("COLL")
160 
161  ABI_CHECK(mkmem/=0,"mkmem==0 not supported anymore!")
162 
163 !Init MPI data
164 ! spaceComm=mpi_enreg%comm_cell
165 ! if(mpi_enreg%paral_kgb==1) spaceComm=mpi_enreg%comm_kpt
166  spaceComm=mpi_enreg%comm_kpt
167  me=mpi_enreg%me_kpt
168 
169 !Check size of cprj
170  my_nspinor=max(1,nspinor/mpi_enreg%nproc_spinor)
171  if (mcprj/=my_nspinor*mband_cprj*mkmem*nsppol) then
172    msg=' wrong size for cprj !'
173    MSG_BUG(msg)
174  end if
175 
176 !Check if cprj is distributed over bands
177  nprocband=(mband/mband_cprj)
178  if (paral_kgb==1.and.nprocband/=mpi_enreg%nproc_band) then
179    msg=' mband/mband_cprj must be equal to nproc_band!'
180    MSG_BUG(msg)
181  end if
182  if (paw_dmft%use_sc_dmft/=0.and.nprocband/=1) then
183    write(msg,'(4a,e14.3,a)') ch10,&
184 &   ' Parallelization over bands is not yet compatible with self-consistency in DMFT !',ch10,&
185 &   ' Calculation is thus restricted to nstep =1.'
186    MSG_WARNING(msg)
187  end if
188 
189  if( usewvl==1 .and. (nprocband/=1)) then
190    write(msg,'(2a)') ch10,&
191 &   '  ERROR: parallelization over bands is not compatible with WAVELETS'
192    MSG_ERROR(msg)
193  end if
194 
195 !Initialise and check dmft variables
196  if(paw_dmft%use_sc_dmft/=0) then
197    nbandc1=(paw_dmft%mbandc-1)*paw_dmft%use_sc_dmft+1
198  else
199    nbandc1=1
200  end if
201 
202 !Size of pawrhoij datastructure
203  nrhoij=size(pawrhoij)
204 
205 !Check if pawrhoij is distributed over atomic sites
206  paral_atom=(nrhoij/=natom.and.mpi_enreg%nproc_atom>1)
207  if (paral_atom.and.nrhoij/=mpi_enreg%my_natom) then
208    msg=' Size of pawrhoij should be natom or my_natom !'
209    MSG_BUG(msg)
210  end if
211 
212 !Allocate temporary cwaveprj storage
213  ABI_DATATYPE_ALLOCATE(cwaveprj,(natom,nspinor))
214  call pawcprj_alloc(cwaveprj,0,dimcprj)
215  if(paw_dmft%use_sc_dmft/=0) then
216    ABI_DATATYPE_ALLOCATE(cwaveprjb,(natom,nspinor))
217    call pawcprj_alloc(cwaveprjb,0,dimcprj)
218  end if
219 
220 !Initialize temporary file (if used)
221  iorder_cprj=0
222 
223 !Build and initialize unpacked rhoij (to be computed here)
224  call pawrhoij_init_unpacked(pawrhoij)
225 
226 !If pawrhoij is MPI-distributed over atomic sites, gather it
227  if (paral_atom) then
228    ABI_DATATYPE_ALLOCATE(pawrhoij_all,(natom))
229  else
230    pawrhoij_all => pawrhoij
231  end if
232 
233 !LOOP OVER SPINS
234  option=1
235  usetimerev=(kptopt>0.and.kptopt<3)
236  bdtot_index=0;ibg=0;jbg=0
237  do isppol=1,nsppol
238 
239 !  LOOP OVER k POINTS
240    do ikpt=1,nkpt
241 
242      nband_k=nband(ikpt+(isppol-1)*nkpt)
243      nband_k_cprj=nband_k/nprocband
244      wtk_k=wtk(ikpt)
245 
246      if(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,isppol,me)) then
247        bdtot_index=bdtot_index+nband_k
248        cycle
249      end if
250 
251      cplex=2;if (istwfk(ikpt)>1) cplex=1
252 
253 !    In case of spinors parallelism, need some extra storage
254      if (mpi_enreg%paral_spinor==1) then
255        nband_k_cprj_used=min(max_nband_cprj,nband_k_cprj)
256        ABI_DATATYPE_ALLOCATE(cprj_tmp,(natom,my_nspinor*nband_k_cprj_used))
257        ABI_DATATYPE_ALLOCATE(cprj_ptr,(natom,   nspinor*nband_k_cprj_used))
258        call pawcprj_alloc(cprj_tmp,0,dimcprj)
259        call pawcprj_alloc(cprj_ptr,0,dimcprj)
260      else
261        cprj_ptr => cprj
262      end if
263 
264 !    LOOP OVER BANDS
265      ib_this_proc=0;jb_this_proc=0
266      do ib=1,nband_k
267        iband=bdtot_index+ib
268 
269 !      Parallelization: treat only some bands
270        if(xmpi_paral==1)then
271          if (paral_kgb==1) then
272            if (mod((ib-1)/mpi_enreg%bandpp,mpi_enreg%nproc_band)/=mpi_enreg%me_band) cycle
273          else
274            if (mpi_enreg%proc_distrb(ikpt,ib,isppol)/=me) cycle
275          end if
276        end if
277        ib_this_proc=ib_this_proc+1
278 
279 !      In case of spinors parallelism, gather cprj because we need both components together
280 !      We do that nband_k_cprj_used by nband_k_cprj_used bands
281        if (mpi_enreg%paral_spinor==1) then
282          jb_this_proc=jb_this_proc+1
283          if (mod(jb_this_proc,nband_k_cprj_used)==1) then
284            ib_this_proc=1
285            nband_k_cprj_read=nband_k_cprj_used
286            if (nband_k_cprj<jb_this_proc+nband_k_cprj_used-1) nband_k_cprj_read=nband_k_cprj-jb_this_proc+1
287            call pawcprj_get(atindx1,cprj_tmp,cprj,natom,jb_this_proc,jbg,ikpt,iorder_cprj,isppol,&
288 &           mband_cprj,mkmem,natom,nband_k_cprj_read,nband_k_cprj,my_nspinor,nsppol,unpaw,&
289 &           mpicomm=mpi_enreg%comm_kpt,proc_distrb=mpi_enreg%proc_distrb)
290            call pawcprj_gather_spin(cprj_tmp,cprj_ptr,natom,nband_k_cprj_read,my_nspinor,nspinor,&
291 &           mpi_enreg%comm_spinor,ierr)
292          end if
293        end if
294 
295 !      DMFT: LOOP ON ADDITIONAL BANDS
296        do ibc1=1,nbandc1
297 !        check if dmft and occupations
298 !        write(std_out,*) 'ib,ibc1          ',ib,ibc1
299 
300 !        DMFT stuff: extract cprj and occupations for additional band
301          if(paw_dmft%use_sc_dmft /= 0) then
302            ib1 = paw_dmft%include_bands(ibc1)
303 !          write(std_out,*) 'use_sc_dmft=1 ib,ib1',ib,ib1
304            iband1 = bdtot_index+ib1
305 !          write(std_out,*) 'ib, ib1          ',paw_dmft%band_in(ib),paw_dmft%band_in(ib1)
306            if(paw_dmft%band_in(ib)) then
307              if(.not.paw_dmft%band_in(ib1))  stop
308              use_nondiag_occup_dmft = 1
309              occup(1) = paw_dmft%occnd(1,ib,ib1,ikpt,isppol)
310              if(nspinor==2) occup(2) = paw_dmft%occnd(2,ib,ib1,ikpt,isppol)
311              if(nspinor==1) occup(2) = zero
312              locc_test = abs(paw_dmft%occnd(1,ib,ib1,ikpt,isppol))+abs(paw_dmft%occnd(2,ib,ib1,ikpt,isppol))>tol8
313 !            write(std_out,*) 'use_sc_dmft=1,band_in(ib)=1, ib,ibc1',ib,ib1,locc_test
314              if (locc_test .or. mkmem == 0) then
315                call pawcprj_get(atindx1,cwaveprjb,cprj_ptr,natom,ib1,ibg,ikpt,iorder_cprj,isppol,&
316 &               mband_cprj,mkmem,natom,1,nband_k_cprj,nspinor,nsppol,unpaw,&
317 &               mpicomm=mpi_enreg%comm_kpt,proc_distrb=mpi_enreg%proc_distrb)
318              end if
319            else
320              use_nondiag_occup_dmft = 0
321              locc_test = (abs(occ(iband))>tol8)
322              occup(1) = occ(iband)
323              if(ibc1 /= 1 .and. .not.(paw_dmft%band_in(ib))) cycle
324            end if
325          else  ! nbandc1=1
326            use_nondiag_occup_dmft=0
327            locc_test = (abs(occ(iband))>tol8)
328            occup(1) = occ(iband)
329          end if
330 
331 !        Extract cprj for current band
332 !        Must read cprj when mkmem=0 (even if unused) to have right pointer inside _PAW file
333          if (locc_test.or.mkmem==0) then
334            call pawcprj_get(atindx1,cwaveprj,cprj_ptr,natom,ib_this_proc,ibg,ikpt,iorder_cprj,isppol,&
335 &           mband_cprj,mkmem,natom,1,nband_k_cprj,nspinor,nsppol,unpaw,&
336 &           mpicomm=mpi_enreg%comm_kpt,proc_distrb=mpi_enreg%proc_distrb)
337          end if
338 
339 !        Accumulate contribution from (occupied) current band
340          if (locc_test) then
341            if(use_nondiag_occup_dmft == 1) then
342              call pawaccrhoij(atindx,cplex,cwaveprj,cwaveprjb,0,isppol,nrhoij,natom,&
343 &             nspinor,occup(1),option,pawrhoij_all,usetimerev,wtk_k,&
344 &             occ_k_2=occup(2))
345            else
346              call pawaccrhoij(atindx,cplex,cwaveprj,cwaveprj ,0,isppol,nrhoij,natom,&
347 &             nspinor,occup(1),option,pawrhoij_all,usetimerev,wtk_k)
348            end if
349          end if
350        end do ! ib1c
351      end do ! ib
352 
353      if (mpi_enreg%paral_spinor==1) then
354        call pawcprj_free(cprj_tmp)
355        call pawcprj_free(cprj_ptr)
356        ABI_DATATYPE_DEALLOCATE(cprj_tmp)
357        ABI_DATATYPE_DEALLOCATE(cprj_ptr)
358      else
359        nullify(cprj_ptr)
360      end if
361 
362      bdtot_index=bdtot_index+nband_k
363      if (mkmem/=0) then
364        if (mpi_enreg%paral_spinor==0) then
365          ibg=ibg+   nspinor*nband_k_cprj
366        else
367          jbg=jbg+my_nspinor*nband_k_cprj
368        end if
369      end if
370 
371    end do ! ikpt
372  end do ! isppol
373 
374 !deallocate temporary cwaveprj/cprj storage
375  call pawcprj_free(cwaveprj)
376  ABI_DATATYPE_DEALLOCATE(cwaveprj)
377  if(paw_dmft%use_sc_dmft/=0) then
378    call pawcprj_free(cwaveprjb)
379    ABI_DATATYPE_DEALLOCATE(cwaveprjb)
380  end if
381 
382 !MPI: need to exchange rhoij_ between procs
383  if (paral_kgb==1.and.nprocband>1) then
384    call pawrhoij_mpisum_unpacked(pawrhoij_all,spaceComm,comm2=mpi_enreg%comm_band)
385  else
386    call pawrhoij_mpisum_unpacked(pawrhoij_all,spaceComm)
387  end if
388 
389 !In case of distribution over atomic sites, dispatch rhoij
390  if (paral_atom) then
391    do iatom=1,nrhoij
392      iatom_tot=mpi_enreg%my_atmtab(iatom)
393      pawrhoij(iatom)%rhoij_(:,:)=pawrhoij_all(iatom_tot)%rhoij_(:,:)
394    end do
395    call pawrhoij_free(pawrhoij_all)
396    ABI_DATATYPE_DEALLOCATE(pawrhoij_all)
397  end if
398 
399 !Print info
400  if (abs(pawprtvol)>=1) then
401    wrt_mode='COLL';if (paral_atom) wrt_mode='PERS'
402    do iatom=1,nrhoij
403      iatom_tot=iatom;if (paral_atom) iatom_tot=mpi_enreg%my_atmtab(iatom)
404      if (pawprtvol>=0.and.iatom_tot/=1.and.iatom_tot/=natom) cycle
405      nsp2=pawrhoij(iatom)%nsppol;if (pawrhoij(iatom)%nspden==4) nsp2=4
406      write(msg, '(4a,i3,a)') ch10," PAW TEST:",ch10,&
407 &     ' ====== Values of RHOIJ in pawmkrhoij (iatom=',iatom_tot,') ======'
408      if (pawrhoij(iatom)%nspden==2.and.pawrhoij(iatom)%nsppol==1) write(msg,'(3a)') trim(msg),ch10,&
409 &     '      (antiferromagnetism case: only one spin component)'
410      call wrtout(std_out,msg,wrt_mode)
411      do isppol=1,nsp2
412        if (pawrhoij(iatom)%nspden/=1) then
413          write(msg,'(3a)') '   Component ',trim(dspin(isppol+2*(pawrhoij(iatom)%nspden/4))),':'
414          call wrtout(std_out,msg,wrt_mode)
415        end if
416        option=2;if (pawrhoij(iatom)%cplex==2.and.pawrhoij(iatom)%nspinor==1) option=1
417        call pawio_print_ij(std_out,pawrhoij(iatom)%rhoij_(:,isppol),pawrhoij(iatom)%lmn2_size,&
418 &       pawrhoij(iatom)%cplex,pawrhoij(iatom)%lmn_size,-1,idum,0,pawprtvol,idum,&
419 &       -1._dp,1,opt_sym=option,mode_paral=wrt_mode)
420      end do
421    end do
422  end if
423 
424  DBG_EXIT("COLL")
425 
426 end subroutine pawmkrhoij