TABLE OF CONTENTS
ABINIT/initrhoij [ Functions ]
NAME
initrhoij
FUNCTION
Initialize PAW rhoij occupancies (in packed storage) from atomic ones
COPYRIGHT
Copyright (C) 1998-2018 ABINIT group (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 . For the initials of contributors, see ~abinit/doc/developers/contributors.txt.
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
53 #if defined HAVE_CONFIG_H 54 #include "config.h" 55 #endif 56 57 #include "abi_common.h" 58 59 subroutine initrhoij(cplex,lexexch,lpawu,my_natom,natom,& 60 & nspden,nspinor,nsppol,ntypat,pawrhoij,pawspnorb,pawtab,spinat,typat,& 61 & ngrhoij,nlmnmix,use_rhoij_,use_rhoijres,& ! optional arguments 62 & mpi_atmtab,comm_atom) ! optional arguments (parallelism) 63 64 use defs_basis 65 use m_profiling_abi 66 use m_errors 67 use m_xmpi, only : xmpi_comm_self 68 69 use m_pawtab, only : pawtab_type 70 use m_pawrhoij, only : pawrhoij_type, pawrhoij_alloc, pawrhoij_get_nspden 71 use m_paral_atom, only : get_my_atmtab, free_my_atmtab 72 73 !This section has been created automatically by the script Abilint (TD). 74 !Do not modify the following lines by hand. 75 #undef ABI_FUNC 76 #define ABI_FUNC 'initrhoij' 77 !End of the abilint section 78 79 implicit none 80 81 !Arguments --------------------------------------------- 82 !scalars 83 integer,intent(in) :: cplex,my_natom,natom,nspden,nspinor,nsppol,ntypat,pawspnorb 84 integer,intent(in),optional :: comm_atom,ngrhoij,nlmnmix,use_rhoij_,use_rhoijres 85 character(len=500) :: message 86 !arrays 87 integer,intent(in) :: lexexch(ntypat),lpawu(ntypat) 88 integer,intent(in) :: typat(natom) 89 integer,optional,target,intent(in) :: mpi_atmtab(:) 90 real(dp),intent(in) :: spinat(3,natom) 91 type(pawrhoij_type),intent(inout) :: pawrhoij(my_natom) 92 type(pawtab_type),intent(in) :: pawtab(ntypat) 93 94 !Local variables --------------------------------------- 95 !Arrays 96 !scalars 97 integer :: iatom,iatom_rhoij,ilmn,ispden,itypat,j0lmn,jl,jlmn,jspden,klmn,klmn1,my_comm_atom 98 integer :: ngrhoij0,nlmnmix0,nselect,nselect1,nspden_rhoij,use_rhoij_0,use_rhoijres0 99 real(dp) :: ratio,ro,roshift,zratio,zz 100 logical :: my_atmtab_allocated,paral_atom,spinat_zero,test_exexch,test_pawu 101 !arrays 102 integer,pointer :: my_atmtab(:) 103 104 !************************************************************************ 105 106 DBG_ENTER("COLL") 107 108 !PAW+U and local exact-exchange restriction 109 do itypat=1,ntypat 110 if (lpawu(itypat)/=lexexch(itypat).and. lpawu(itypat)/=-1.and.lexexch(itypat)/=-1) then 111 message = ' lpawu must be equal to lexexch !' 112 MSG_ERROR(message) 113 end if 114 end do 115 116 !Set up parallelism over atoms 117 paral_atom=(present(comm_atom).and.(my_natom/=natom)) 118 nullify(my_atmtab);if (present(mpi_atmtab)) my_atmtab => mpi_atmtab 119 my_comm_atom=xmpi_comm_self;if (present(comm_atom)) my_comm_atom=comm_atom 120 call get_my_atmtab(my_comm_atom,my_atmtab,my_atmtab_allocated,paral_atom,natom,my_natom_ref=my_natom) 121 122 nspden_rhoij=pawrhoij_get_nspden(nspden,nspinor,pawspnorb) 123 ratio=one;if (nspden_rhoij==2) ratio=half 124 spinat_zero=all(abs(spinat(:,:))<tol10) 125 126 if (my_natom>0) then 127 ngrhoij0=0;if (present(ngrhoij)) ngrhoij0=ngrhoij 128 nlmnmix0=0;if (present(nlmnmix)) nlmnmix0=nlmnmix 129 use_rhoij_0=0;if (present(use_rhoij_)) use_rhoij_0=use_rhoij_ 130 use_rhoijres0=0;if (present(use_rhoijres)) use_rhoijres0=use_rhoijres 131 if (paral_atom) then 132 call pawrhoij_alloc(pawrhoij,cplex,nspden_rhoij,nspinor,nsppol,typat,& 133 & ngrhoij=ngrhoij0,nlmnmix=nlmnmix0,use_rhoij_=use_rhoij_0,use_rhoijres=use_rhoijres0,& 134 & pawtab=pawtab,comm_atom=my_comm_atom,mpi_atmtab=my_atmtab) 135 else 136 call pawrhoij_alloc(pawrhoij,cplex,nspden_rhoij,nspinor,nsppol,typat,pawtab=pawtab,& 137 & ngrhoij=ngrhoij0,nlmnmix=nlmnmix0,use_rhoij_=use_rhoij_0,use_rhoijres=use_rhoijres0) 138 end if 139 end if 140 141 do iatom_rhoij=1,my_natom 142 iatom=iatom_rhoij;if (paral_atom) iatom=my_atmtab(iatom_rhoij) 143 itypat=typat(iatom) 144 nselect=0 145 146 ! Determine Z (trace of rhoij0 or part of it) 147 zz=zero 148 do jlmn=1,pawtab(itypat)%lmn_size 149 jl=pawtab(itypat)%indlmn(1,jlmn) 150 j0lmn=jlmn*(jlmn-1)/2 151 test_pawu=(lpawu(itypat)==-1.or.lpawu(itypat)==jl) 152 test_exexch=(lexexch(itypat)==-1.or.lexexch(itypat)==jl) 153 do ilmn=1,jlmn 154 klmn=j0lmn+ilmn 155 if ((ilmn==jlmn).and.test_pawu.and.test_exexch) & 156 & zz=zz+pawtab(itypat)%rhoij0(klmn) 157 end do 158 end do 159 160 ! Compute rhoij from tabulated value and magnetization 161 do ispden=1,nspden_rhoij 162 163 zratio=zero 164 roshift=one 165 ratio=one 166 if (nspden_rhoij==2) then 167 ratio=half 168 if ((spinat(3,iatom)>zero.and.ispden==1).or.& 169 & (spinat(3,iatom)<zero.and.ispden==2)) then 170 if(abs(zz)>tol12)then 171 zratio=two*abs(spinat(3,iatom))/zz 172 else 173 zratio=zero 174 end if 175 end if 176 else if (nspden_rhoij==4.and.ispden>=2) then 177 roshift=zero 178 if(abs(zz)>tol12)then 179 zratio=spinat(ispden-1,iatom)/zz 180 else 181 zratio=zero 182 end if 183 end if 184 185 nselect=0;nselect1=1-cplex 186 do jlmn=1,pawtab(itypat)%lmn_size 187 jl=pawtab(itypat)%indlmn(1,jlmn) 188 j0lmn=jlmn*(jlmn-1)/2 189 test_pawu=(lpawu(itypat)==-1.or.lpawu(itypat)==jl) 190 test_exexch=(lexexch(itypat)==-1.or.lexexch(itypat)==jl) 191 do ilmn=1,jlmn 192 klmn=j0lmn+ilmn 193 ro=pawtab(itypat)%rhoij0(klmn) 194 if ((ilmn==jlmn).and.test_pawu.and.test_exexch) then 195 ro=ro*ratio*(roshift+zratio) 196 else 197 ro=ro*ratio*roshift 198 end if 199 200 klmn1=cplex*(klmn-1)+1 201 if (abs(ro)>tol10) then 202 pawrhoij(iatom_rhoij)%rhoijp(klmn1,ispden)=ro 203 else 204 pawrhoij(iatom_rhoij)%rhoijp(klmn1,ispden)=zero 205 end if 206 207 if (ispden==nspden_rhoij) then 208 if (any(abs(pawrhoij(iatom_rhoij)%rhoijp(klmn1,:))>tol10)) then 209 nselect=nselect+1;nselect1=nselect1+cplex 210 pawrhoij(iatom_rhoij)%rhoijselect(nselect)=klmn 211 do jspden=1,nspden_rhoij 212 pawrhoij(iatom_rhoij)%rhoijp(nselect1,jspden)=pawrhoij(iatom_rhoij)%rhoijp(klmn1,jspden) 213 end do 214 end if 215 end if 216 217 end do 218 end do 219 220 end do 221 pawrhoij(iatom_rhoij)%nrhoijsel=nselect 222 223 ! Non-collinear magnetism: avoid zero magnetization, because it produces numerical instabilities 224 ! Add a small real to the magnetization ; not yet activated => must be tested. 225 ! if (pawrhoij(iatom_rhoij)%nspden==4.and.spinat_zero) then 226 ! pawrhoij(iatom_rhoij)%rhoijp(:,4)=pawrhoij(iatom_rhoij)%rhoijp(:,4)+tol10 227 ! end if 228 229 end do ! iatom_rhoij 230 231 !Destroy atom table used for parallelism 232 call free_my_atmtab(my_atmtab,my_atmtab_allocated) 233 234 DBG_EXIT("COLL") 235 236 end subroutine initrhoij