TABLE OF CONTENTS
- m_paw_occupancies/initrhoij
- m_paw_occupancies/m_paw_occupancies
- m_paw_occupancies/pawaccrhoij
- m_paw_occupancies/pawmkrhoij
m_paw_occupancies/initrhoij [ Functions ]
[ Top ] [ m_paw_occupancies ] [ Functions ]
NAME
initrhoij
FUNCTION
Initialize PAW rhoij occupancies (in packed storage) from atomic ones
INPUTS
cpxocc=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) qphase=2 if rhoij have a exp(iqR) phase, 1 if not (typical use: response function at q<>0) 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
SOURCE
1173 subroutine initrhoij(cpxocc,lexexch,lpawu,my_natom,natom,nspden,nspinor,nsppol,& 1174 & ntypat,pawrhoij,pawspnorb,pawtab,qphase,spinat,typat,& 1175 & ngrhoij,nlmnmix,use_rhoij_,use_rhoijres,& ! optional arguments 1176 & mpi_atmtab,comm_atom) ! optional arguments (parallelism) 1177 1178 !Arguments --------------------------------------------- 1179 !scalars 1180 integer,intent(in) :: cpxocc,my_natom,natom,nspden,nspinor,nsppol,ntypat,pawspnorb,qphase 1181 integer,intent(in),optional :: comm_atom,ngrhoij,nlmnmix,use_rhoij_,use_rhoijres 1182 character(len=500) :: message 1183 !arrays 1184 integer,intent(in) :: lexexch(ntypat),lpawu(ntypat) 1185 integer,intent(in) :: typat(natom) 1186 integer,optional,target,intent(in) :: mpi_atmtab(:) 1187 real(dp),intent(in) :: spinat(3,natom) 1188 type(pawrhoij_type),intent(inout) :: pawrhoij(my_natom) 1189 type(pawtab_type),intent(in) :: pawtab(ntypat) 1190 1191 !Local variables --------------------------------------- 1192 !Arrays 1193 !scalars 1194 integer :: cplex_rhoij,iatom,iatom_rhoij,ilmn,ispden,itypat,j0lmn,jl,jlmn,jspden 1195 integer :: klmn,klmn1,ln,lnspinat0,my_comm_atom 1196 integer :: ngrhoij0,nlmnmix0,nselect,nselect1,nspden_rhoij,qphase_rhoij 1197 integer :: use_rhoij_0,use_rhoijres0 1198 real(dp) :: ratio,ro,roshift,zratio,zz 1199 logical :: my_atmtab_allocated,paral_atom,spinat_zero,test_exexch,test_pawu,test_lnspinat 1200 !arrays 1201 integer,pointer :: my_atmtab(:),lnspinat(:) 1202 real(dp),allocatable :: occ(:) 1203 1204 !************************************************************************ 1205 1206 DBG_ENTER("COLL") 1207 1208 !PAW+U and local exact-exchange restriction 1209 do itypat=1,ntypat 1210 if (lpawu(itypat)/=lexexch(itypat).and. lpawu(itypat)/=-1.and.lexexch(itypat)/=-1) then 1211 message = ' lpawu must be equal to lexexch !' 1212 ABI_ERROR(message) 1213 end if 1214 end do 1215 1216 !Set up parallelism over atoms 1217 paral_atom=(present(comm_atom).and.(my_natom/=natom)) 1218 nullify(my_atmtab);if (present(mpi_atmtab)) my_atmtab => mpi_atmtab 1219 my_comm_atom=xmpi_comm_self;if (present(comm_atom)) my_comm_atom=comm_atom 1220 call get_my_atmtab(my_comm_atom,my_atmtab,my_atmtab_allocated,paral_atom,natom,my_natom_ref=my_natom) 1221 1222 call pawrhoij_inquire_dim(cplex_rhoij=cplex_rhoij,qphase_rhoij=qphase_rhoij,nspden_rhoij=nspden_rhoij,& 1223 & nspden=nspden,spnorb=pawspnorb,cpxocc=cpxocc,cplex=qphase) 1224 1225 ratio=one;if (nspden_rhoij==2) ratio=half 1226 spinat_zero=(all(abs(spinat(:,:))<tol10).or.(nspden_rhoij==4.and.nspden==1)) 1227 1228 ! if (nspden_rhoij==4.and.nspden==1.and.(.not.spinat_zero)) then 1229 ! write(message,'(5a)') 'You are performing a unpolarized calculation (nspden==1)',ch10,& 1230 !& 'but you start with a magnetization on atom (spinat/=0).',ch10,& 1231 !& 'This is not expected and my produce unphysical results!' 1232 ! ABI_WARNING(message) 1233 ! end if 1234 1235 if (my_natom>0) then 1236 ngrhoij0=0;if (present(ngrhoij)) ngrhoij0=ngrhoij 1237 nlmnmix0=0;if (present(nlmnmix)) nlmnmix0=nlmnmix 1238 use_rhoij_0=0;if (present(use_rhoij_)) use_rhoij_0=use_rhoij_ 1239 use_rhoijres0=0;if (present(use_rhoijres)) use_rhoijres0=use_rhoijres 1240 if (paral_atom) then 1241 call pawrhoij_alloc(pawrhoij,cplex_rhoij,nspden_rhoij,nspinor,nsppol,typat,& 1242 & ngrhoij=ngrhoij0,nlmnmix=nlmnmix0,use_rhoij_=use_rhoij_0,use_rhoijres=use_rhoijres0,& 1243 & qphase=qphase_rhoij,pawtab=pawtab,comm_atom=my_comm_atom,mpi_atmtab=my_atmtab) 1244 else 1245 call pawrhoij_alloc(pawrhoij,cplex_rhoij,nspden_rhoij,nspinor,nsppol,typat,qphase=qphase_rhoij,& 1246 & pawtab=pawtab,ngrhoij=ngrhoij0,nlmnmix=nlmnmix0,use_rhoij_=use_rhoij_0,use_rhoijres=use_rhoijres0) 1247 end if 1248 end if 1249 1250 do iatom_rhoij=1,my_natom 1251 iatom=iatom_rhoij;if (paral_atom) iatom=my_atmtab(iatom_rhoij) 1252 itypat=typat(iatom) 1253 nselect=0 1254 ABI_MALLOC(lnspinat,(pawtab(itypat)%basis_size)) 1255 lnspinat=-1 1256 ! Determine occupancies of each orbital 1257 if (nspden_rhoij==2) then 1258 ABI_MALLOC(occ,(pawtab(itypat)%basis_size)) 1259 occ=zero 1260 do jlmn=1,pawtab(itypat)%lmn_size 1261 ln=pawtab(itypat)%indlmn(5,jlmn) 1262 klmn=jlmn*(jlmn+1)/2 1263 occ(ln)=occ(ln)+pawtab(itypat)%rhoij0(klmn) 1264 end do 1265 do ln=1,pawtab(itypat)%basis_size 1266 if(pawtab(itypat)%orbitals(ln)==0.and.occ(ln)==1) lnspinat(ln)=ln 1267 if(pawtab(itypat)%orbitals(ln)==1.and.(occ(ln)>=1.and.occ(ln)<=5)) lnspinat(ln)=ln 1268 if(pawtab(itypat)%orbitals(ln)==2.and.(occ(ln)>=1.and.occ(ln)<=9)) lnspinat(ln)=ln 1269 if(pawtab(itypat)%orbitals(ln)==3.and.(occ(ln)>=1.and.occ(ln)<=13)) lnspinat(ln)=ln 1270 end do 1271 ABI_FREE(occ) 1272 end if 1273 lnspinat0=maxval(lnspinat) 1274 lnspinat0=-1 1275 1276 ! Determine Z (trace of rhoij0 or part of it) 1277 zz=zero 1278 do jlmn=1,pawtab(itypat)%lmn_size 1279 jl=pawtab(itypat)%indlmn(1,jlmn) 1280 ln=pawtab(itypat)%indlmn(5,jlmn) 1281 j0lmn=jlmn*(jlmn-1)/2 1282 test_lnspinat=(lnspinat0==-1.or.lnspinat(ln)==ln) 1283 test_pawu=(lpawu(itypat)==-1.or.lpawu(itypat)==jl) 1284 test_exexch=(lexexch(itypat)==-1.or.lexexch(itypat)==jl) 1285 do ilmn=1,jlmn 1286 klmn=j0lmn+ilmn 1287 if ((ilmn==jlmn).and.test_pawu.and.test_exexch.and.test_lnspinat) & 1288 & zz=zz+pawtab(itypat)%rhoij0(klmn) 1289 end do 1290 end do 1291 1292 ! Compute rhoij from tabulated value and magnetization 1293 do ispden=1,nspden_rhoij 1294 1295 zratio=zero 1296 roshift=one 1297 ratio=one 1298 if (nspden_rhoij==2) then 1299 ratio=half 1300 if ((spinat(3,iatom)>tol12.and.ispden==1).or.& 1301 & (spinat(3,iatom)<tol12.and.ispden==2)) then 1302 if(abs(zz)>tol12)then 1303 zratio=two*abs(spinat(3,iatom))/zz 1304 else 1305 zratio=zero 1306 end if 1307 end if 1308 else if (nspden_rhoij==4.and.ispden>=2) then 1309 roshift=zero 1310 if(abs(zz)>tol12.and.(.not.spinat_zero)) then 1311 zratio=spinat(ispden-1,iatom)/zz 1312 else 1313 zratio=zero 1314 end if 1315 end if 1316 1317 nselect=0;nselect1=1-cpxocc 1318 do jlmn=1,pawtab(itypat)%lmn_size 1319 jl=pawtab(itypat)%indlmn(1,jlmn) 1320 ln=pawtab(itypat)%indlmn(5,jlmn) 1321 j0lmn=jlmn*(jlmn-1)/2 1322 test_lnspinat=(lnspinat0==-1.or.lnspinat(ln)==ln) 1323 test_pawu=(lpawu(itypat)==-1.or.lpawu(itypat)==jl) 1324 test_exexch=(lexexch(itypat)==-1.or.lexexch(itypat)==jl) 1325 do ilmn=1,jlmn 1326 klmn=j0lmn+ilmn 1327 ro=pawtab(itypat)%rhoij0(klmn) 1328 if ((ilmn==jlmn).and.test_pawu.and.test_exexch.and.test_lnspinat) then 1329 ro=ro*ratio*(roshift+zratio) 1330 else 1331 ro=ro*ratio*roshift 1332 end if 1333 1334 klmn1=cpxocc*(klmn-1)+1 1335 if (abs(ro)>tol10) then 1336 pawrhoij(iatom_rhoij)%rhoijp(klmn1,ispden)=ro 1337 else 1338 pawrhoij(iatom_rhoij)%rhoijp(klmn1,ispden)=zero 1339 end if 1340 1341 if (ispden==nspden_rhoij) then 1342 if (any(abs(pawrhoij(iatom_rhoij)%rhoijp(klmn1,:))>tol10)) then 1343 nselect=nselect+1;nselect1=nselect1+cpxocc 1344 pawrhoij(iatom_rhoij)%rhoijselect(nselect)=klmn 1345 do jspden=1,nspden_rhoij 1346 pawrhoij(iatom_rhoij)%rhoijp(nselect1,jspden)=pawrhoij(iatom_rhoij)%rhoijp(klmn1,jspden) 1347 end do 1348 end if 1349 end if 1350 1351 end do 1352 end do 1353 1354 end do 1355 pawrhoij(iatom_rhoij)%nrhoijsel=nselect 1356 if (nselect<pawrhoij(iatom_rhoij)%lmn2_size) & 1357 & pawrhoij(iatom_rhoij)%rhoijselect(nselect+1:pawrhoij(iatom_rhoij)%lmn2_size)=0 1358 1359 ! Non-collinear magnetism: avoid zero magnetization, because it produces numerical instabilities 1360 ! Add a small real to the magnetization ; not yet activated => must be tested. 1361 ! if (pawrhoij(iatom_rhoij)%nspden==4.and.spinat_zero) then 1362 ! pawrhoij(iatom_rhoij)%rhoijp(:,4)=pawrhoij(iatom_rhoij)%rhoijp(:,4)+tol10 1363 ! end if 1364 ABI_FREE(lnspinat) 1365 end do ! iatom_rhoij 1366 1367 !Destroy atom table used for parallelism 1368 call free_my_atmtab(my_atmtab,my_atmtab_allocated) 1369 1370 DBG_EXIT("COLL") 1371 1372 end subroutine initrhoij
m_paw_occupancies/m_paw_occupancies [ Modules ]
NAME
m_paw_occupancies
FUNCTION
This module contains routines related to the computation of PAW on-site occupancies (rhoij).
COPYRIGHT
Copyright (C) 2018-2024 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
16 #if defined HAVE_CONFIG_H 17 #include "config.h" 18 #endif 19 20 #include "abi_common.h" 21 22 MODULE m_paw_occupancies 23 24 use defs_basis 25 use m_abicore 26 use m_errors 27 use m_xmpi 28 29 use defs_abitypes, only : MPI_type 30 use m_pawtab, only : pawtab_type 31 use m_pawrhoij, only : pawrhoij_type,pawrhoij_init_unpacked,pawrhoij_mpisum_unpacked, & 32 & pawrhoij_alloc,pawrhoij_free,pawrhoij_inquire_dim 33 use m_pawcprj, only : pawcprj_type,pawcprj_alloc,pawcprj_get, & 34 & pawcprj_gather_spin, pawcprj_free, pawcprj_mpi_send, & 35 & pawcprj_mpi_recv, pawcprj_copy, pawcprj_unpack, pawcprj_pack 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
use_timerev=.TRUE. if time-reversal symmetry is used (WF(-k)=Conjg[WF(k)])
use_zeromag=.TRUE. if rhoij "magnetization" is enforced to be zero
Applies only when nspden_rhoij=4 (note: only the real part is set to zero)
wtk_k=weight assigned to current k-point
[occ_k_2]=??
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)}
SOURCE
566 subroutine pawaccrhoij(atindx,cplex,cwaveprj,cwaveprj1,ipert,isppol,my_natom,natom,& 567 & nspinor,occ_k,option,pawrhoij,use_timerev,use_zeromag,wtk_k,occ_k_2, & 568 & comm_atom,mpi_atmtab ) ! optional (parallelism) 569 570 !Arguments --------------------------------------------- 571 !scalars 572 integer,intent(in) :: cplex,ipert,isppol,my_natom,natom,nspinor,option 573 integer,optional,intent(in) :: comm_atom 574 logical,intent(in) :: use_timerev,use_zeromag 575 real(dp),intent(in) :: occ_k,wtk_k 576 real(dp),optional,intent(in) :: occ_k_2 577 !arrays 578 integer,intent(in) :: atindx(natom) 579 integer,optional,target,intent(in) :: mpi_atmtab(:) 580 type(pawcprj_type),intent(in) :: cwaveprj(natom,nspinor),cwaveprj1(natom,nspinor) 581 type(pawrhoij_type),intent(inout) :: pawrhoij(my_natom) 582 583 !Local variables --------------------------------------- 584 !scalars 585 integer :: cplex_rhoij,iatm,iatom,iatom1,ilmn,iplex,iq0,j0lmn,jlmn,klmn,klmn_im,klmn_re 586 integer :: mu,my_comm_atom,ncpgr,nspden_rhoij 587 logical :: compute_impart,compute_impart_cplex,substract_diagonal 588 logical :: my_atmtab_allocated,paral_atom 589 real(dp) :: ro11_im,ro11_re,ro12_im,ro12_re,ro21_im,ro21_re,ro22_im,ro22_re,weight,weight_2 590 character(len=500) :: message 591 !arrays 592 integer,pointer :: my_atmtab(:) 593 real(dp) :: cpi0(2,nspinor),cpi1(2,nspinor),cpj0(2,nspinor),cpj1(2,nspinor) 594 real(dp) :: dcpi0(2,nspinor,9),dcpj0(2,nspinor,9) 595 596 ! *********************************************************************** 597 598 DBG_ENTER("COLL") 599 600 if (my_natom==0) return 601 602 ncpgr=0 603 if (option==2.and.(ipert<=natom.or.ipert==natom+3.or.ipert==natom+4)) ncpgr=1 604 if (option==3) ncpgr=cwaveprj(1,1)%ncpgr 605 606 !Tests 607 if(option==2.and.(ipert==natom+1.or.ipert==natom+10.or.ipert==natom+11)) then 608 message = 'Not relevant for ipert=natom+1 or ipert=natom+10 or ipert=natom+11!' 609 ABI_BUG(message) 610 end if 611 if(option==2.and.cwaveprj(1,1)%ncpgr<ncpgr) then 612 message = 'Error on cwaveprj1 factors derivatives!' 613 ABI_BUG(message) 614 end if 615 if(option==3.and.cwaveprj(1,1)%ncpgr/=ncpgr) then 616 message = 'Error on cwaveprj factors derivatives!' 617 ABI_BUG(message) 618 end if 619 if (pawrhoij(1)%qphase==2.and.option/=2) then 620 message = 'pawaccrhoij: qphase=2 only allowed with option=2 (1st-order rhoij)!' 621 ABI_BUG(message) 622 end if 623 624 !Set up parallelism over atoms 625 paral_atom=(present(comm_atom).and.(my_natom/=natom)) 626 nullify(my_atmtab);if (present(mpi_atmtab)) my_atmtab => mpi_atmtab 627 my_comm_atom=xmpi_comm_self;if (present(comm_atom)) my_comm_atom=comm_atom 628 call get_my_atmtab(my_comm_atom,my_atmtab,my_atmtab_allocated,paral_atom,natom,& 629 & my_natom_ref=my_natom) 630 631 weight=wtk_k*occ_k 632 weight_2=zero 633 if(present(occ_k_2)) then 634 if (cplex==1) then 635 !DMFT need complex cprj 636 message='DMFT computation must be done with complex cprj. Check that istwfk = 1!' 637 ABI_ERROR(message) 638 end if 639 weight_2=wtk_k*occ_k_2 640 end if 641 if (pawrhoij(1)%nspden==2.and.pawrhoij(1)%nsppol==1.and.nspinor==1) weight=half*weight 642 if (pawrhoij(1)%nspden==2.and.pawrhoij(1)%nsppol==1.and.nspinor==1.and.present(occ_k_2)) weight_2=half*weight_2 643 644 if (option==1) then 645 646 ! ================================================================== 647 ! === OPTION 1: Accumulate (n,k) contribution to rhoij ============= 648 ! ================================================================== 649 compute_impart=((.not.use_timerev).and.(pawrhoij(1)%cplex_rhoij==2)) 650 compute_impart_cplex=((compute_impart).and.(cplex==2)) 651 if (nspinor==1) then 652 cplex_rhoij=pawrhoij(1)%cplex_rhoij 653 if (cplex_rhoij==1) then 654 do iatom=1,my_natom 655 iatom1=iatom;if (paral_atom) iatom1=my_atmtab(iatom) 656 iatm=atindx(iatom1) 657 do jlmn=1,pawrhoij(iatom)%lmn_size 658 j0lmn=jlmn*(jlmn-1)/2 659 cpj0(1:cplex,1)=cwaveprj(iatm,1)%cp(1:cplex,jlmn) 660 do ilmn=1,jlmn 661 klmn=j0lmn+ilmn 662 cpi0(1:cplex,1)=cwaveprj1(iatm,1)%cp(1:cplex,ilmn) 663 ro11_re=zero 664 do iplex=1,cplex 665 ro11_re=ro11_re+cpi0(iplex,1)*cpj0(iplex,1) 666 end do 667 pawrhoij(iatom)%rhoij_(klmn,isppol)=pawrhoij(iatom)%rhoij_(klmn,isppol)+weight*ro11_re 668 if (present(occ_k_2)) then 669 ro11_im=cpi0(1,1)*cpj0(2,1)-cpi0(2,1)*cpj0(1,1) 670 pawrhoij(iatom)%rhoij_(klmn,isppol)=pawrhoij(iatom)%rhoij_(klmn,isppol)-weight_2*ro11_im 671 end if 672 end do 673 end do 674 end do 675 else 676 do iatom=1,my_natom 677 iatom1=iatom;if (paral_atom) iatom1=my_atmtab(iatom) 678 iatm=atindx(iatom1) 679 do jlmn=1,pawrhoij(iatom)%lmn_size 680 j0lmn=jlmn*(jlmn-1)/2 681 cpj0(1:cplex,1)=cwaveprj(iatm,1)%cp(1:cplex,jlmn) 682 do ilmn=1,jlmn 683 klmn=j0lmn+ilmn 684 klmn_re=cplex_rhoij*(klmn-1)+1 685 cpi0(1:cplex,1)=cwaveprj1(iatm,1)%cp(1:cplex,ilmn) 686 ro11_re=zero 687 do iplex=1,cplex 688 ro11_re=ro11_re+cpi0(iplex,1)*cpj0(iplex,1) 689 end do 690 pawrhoij(iatom)%rhoij_(klmn_re,isppol)=pawrhoij(iatom)%rhoij_(klmn_re,isppol)+weight*ro11_re 691 if (present(occ_k_2)) then 692 ro11_im=cpi0(1,1)*cpj0(2,1)-cpi0(2,1)*cpj0(1,1) 693 pawrhoij(iatom)%rhoij_(klmn_re,isppol)=pawrhoij(iatom)%rhoij_(klmn_re,isppol)-weight_2*ro11_im 694 end if 695 if (compute_impart_cplex) then 696 klmn_im=klmn_re+1 697 ro11_im=cpi0(1,1)*cpj0(2,1)-cpi0(2,1)*cpj0(1,1) 698 pawrhoij(iatom)%rhoij_(klmn_im,isppol)=pawrhoij(iatom)%rhoij_(klmn_im,isppol)+weight*ro11_im 699 if (present(occ_k_2)) then 700 pawrhoij(iatom)%rhoij_(klmn_im,isppol)=pawrhoij(iatom)%rhoij_(klmn_im,isppol)+weight_2*ro11_re 701 end if 702 end if 703 end do 704 end do 705 end do 706 end if 707 else ! nspinor=2 708 do iatom=1,my_natom 709 iatom1=iatom;if (paral_atom) iatom1=my_atmtab(iatom) 710 iatm=atindx(iatom1) 711 cplex_rhoij=pawrhoij(iatom)%cplex_rhoij 712 nspden_rhoij=pawrhoij(iatom)%nspden 713 do jlmn=1,pawrhoij(iatom)%lmn_size 714 j0lmn=jlmn*(jlmn-1)/2 715 cpj0(1:cplex,1)=cwaveprj(iatm,1)%cp(1:cplex,jlmn) 716 cpj0(1:cplex,2)=cwaveprj(iatm,2)%cp(1:cplex,jlmn) 717 do ilmn=1,jlmn 718 klmn=j0lmn+ilmn 719 klmn_re=cplex_rhoij*(klmn-1)+1 720 cpi0(1:cplex,1)=cwaveprj1(iatm,1)%cp(1:cplex,ilmn) 721 cpi0(1:cplex,2)=cwaveprj1(iatm,2)%cp(1:cplex,ilmn) 722 ro11_re=zero;ro22_re=zero 723 ro12_re=zero;ro21_re=zero 724 ro12_im=zero;ro21_im=zero 725 do iplex=1,cplex 726 ro11_re=ro11_re+cpi0(iplex,1)*cpj0(iplex,1) 727 ro22_re=ro22_re+cpi0(iplex,2)*cpj0(iplex,2) 728 end do 729 pawrhoij(iatom)%rhoij_(klmn_re,1)=pawrhoij(iatom)%rhoij_(klmn_re,1)+weight*(ro11_re+ro22_re) 730 if (nspden_rhoij>1) then 731 do iplex=1,cplex 732 ro12_re=ro12_re+cpi0(iplex,2)*cpj0(iplex,1) 733 ro21_re=ro21_re+cpi0(iplex,1)*cpj0(iplex,2) 734 end do 735 if (.not.use_zeromag) then 736 pawrhoij(iatom)%rhoij_(klmn_re,4)=pawrhoij(iatom)%rhoij_(klmn_re,4)+weight*(ro11_re-ro22_re) 737 pawrhoij(iatom)%rhoij_(klmn_re,2)=pawrhoij(iatom)%rhoij_(klmn_re,2)+weight*(ro12_re+ro21_re) 738 end if 739 if (cplex==2) then 740 !Important note: the present implementation follows eq(15) in Hobbs et al, PRB 62, 11556(2000) 741 ! rho^alpha,beta_ij = Sum[<Psi^beta|pi><pj|Psi^alpha] (alpha and beta exponents inverted) 742 ro12_im=cpi0(1,2)*cpj0(2,1)-cpi0(2,2)*cpj0(1,1) 743 ro21_im=cpi0(1,1)*cpj0(2,2)-cpi0(2,1)*cpj0(1,2) 744 if (.not.use_zeromag) then 745 pawrhoij(iatom)%rhoij_(klmn_re,3)=pawrhoij(iatom)%rhoij_(klmn_re,3)+weight*(ro21_im-ro12_im) 746 end if 747 end if 748 end if 749 if (present(occ_k_2)) then 750 ro11_im=cpi0(1,1)*cpj0(2,1)-cpi0(2,1)*cpj0(1,1) 751 ro22_im=cpi0(1,2)*cpj0(2,2)-cpi0(2,2)*cpj0(1,2) 752 pawrhoij(iatom)%rhoij_(klmn_re,1)=pawrhoij(iatom)%rhoij_(klmn_re,1)+weight_2*(-ro11_im-ro22_im) 753 pawrhoij(iatom)%rhoij_(klmn_re,2)=pawrhoij(iatom)%rhoij_(klmn_re,2)+weight_2*(-ro21_im-ro12_im) 754 pawrhoij(iatom)%rhoij_(klmn_re,3)=pawrhoij(iatom)%rhoij_(klmn_re,3)+weight_2*(-ro12_re+ro21_re) 755 pawrhoij(iatom)%rhoij_(klmn_re,4)=pawrhoij(iatom)%rhoij_(klmn_re,4)+weight_2*(-ro11_im+ro22_im) 756 end if 757 if (compute_impart) then 758 klmn_im=klmn_re+1 759 if (nspden_rhoij>1) pawrhoij(iatom)%rhoij_(klmn_im,3)=pawrhoij(iatom)%rhoij_(klmn_im,3)+weight*(ro12_re-ro21_re) 760 if (cplex==2) then 761 ro11_im=cpi0(1,1)*cpj0(2,1)-cpi0(2,1)*cpj0(1,1) 762 ro22_im=cpi0(1,2)*cpj0(2,2)-cpi0(2,2)*cpj0(1,2) 763 pawrhoij(iatom)%rhoij_(klmn_im,1)=pawrhoij(iatom)%rhoij_(klmn_im,1)+weight*(ro11_im+ro22_im) 764 if (nspden_rhoij>1) then 765 pawrhoij(iatom)%rhoij_(klmn_im,4)=pawrhoij(iatom)%rhoij_(klmn_im,4)+weight*(ro11_im-ro22_im) 766 pawrhoij(iatom)%rhoij_(klmn_im,2)=pawrhoij(iatom)%rhoij_(klmn_im,2)+weight*(ro12_im+ro21_im) 767 end if 768 if (present(occ_k_2)) then 769 pawrhoij(iatom)%rhoij_(klmn_im,1)=pawrhoij(iatom)%rhoij_(klmn_im,1)+weight_2*( ro11_re+ro22_re) 770 pawrhoij(iatom)%rhoij_(klmn_im,2)=pawrhoij(iatom)%rhoij_(klmn_im,2)+weight_2*( ro21_re+ro12_re) 771 pawrhoij(iatom)%rhoij_(klmn_im,3)=pawrhoij(iatom)%rhoij_(klmn_im,3)+weight_2*(-ro12_im+ro21_im) 772 pawrhoij(iatom)%rhoij_(klmn_im,4)=pawrhoij(iatom)%rhoij_(klmn_im,4)+weight_2*( ro11_re-ro22_re) 773 end if 774 end if 775 end if 776 end do 777 end do 778 end do 779 end if 780 781 else if (option==2) then 782 783 ! ================================================================== 784 ! === OPTION 2: Accumulate (n,k) contribution to 1st-order rhoij === 785 ! ================================================================== 786 787 ! Accumulate (n,k) contribution to rhoij1 788 ! due to derivative of wave-function 789 compute_impart=(pawrhoij(1)%qphase==2) 790 compute_impart_cplex=(compute_impart.and.(cplex==2)) 791 if (nspinor==1) then 792 do iatom=1,my_natom 793 iatom1=iatom;if (paral_atom) iatom1=my_atmtab(iatom) 794 iatm=atindx(iatom1) 795 cplex_rhoij=pawrhoij(iatom)%cplex_rhoij 796 iq0=cplex_rhoij*pawrhoij(iatom)%lmn2_size 797 do jlmn=1,pawrhoij(iatom)%lmn_size 798 j0lmn=jlmn*(jlmn-1)/2 799 cpj0(1:2,1)=cwaveprj (iatm,1)%cp(1:2,jlmn) 800 cpj1(1:2,1)=cwaveprj1(iatm,1)%cp(1:2,jlmn) 801 do ilmn=1,jlmn 802 klmn=j0lmn+ilmn 803 klmn_re=cplex_rhoij*(klmn-1)+1 804 cpi0(1:2,1)=cwaveprj (iatm,1)%cp(1:2,ilmn) 805 cpi1(1:2,1)=cwaveprj1(iatm,1)%cp(1:2,ilmn) 806 ro11_re=zero 807 do iplex=1,cplex 808 ro11_re=ro11_re+cpi0(iplex,1)*cpj1(iplex,1)+cpj0(iplex,1)*cpi1(iplex,1) 809 end do 810 pawrhoij(iatom)%rhoij_(klmn_re,isppol)=pawrhoij(iatom)%rhoij_(klmn_re,isppol)+weight*ro11_re 811 if (compute_impart_cplex) then 812 klmn_im=klmn_re+iq0 813 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) 814 pawrhoij(iatom)%rhoij_(klmn_im,isppol)=pawrhoij(iatom)%rhoij_(klmn_im,isppol)+weight*ro11_im 815 end if 816 end do 817 end do 818 end do 819 else ! nspinor=2 820 do iatom=1,my_natom 821 iatom1=iatom;if (paral_atom) iatom1=my_atmtab(iatom) 822 iatm=atindx(iatom1) 823 cplex_rhoij=pawrhoij(iatom)%cplex_rhoij 824 nspden_rhoij=pawrhoij(iatom)%nspden 825 iq0=cplex_rhoij*pawrhoij(iatom)%lmn2_size 826 do jlmn=1,pawrhoij(iatom)%lmn_size 827 j0lmn=jlmn*(jlmn-1)/2 828 cpj0(1:2,1)=cwaveprj (iatm,1)%cp(1:2,jlmn) 829 cpj0(1:2,2)=cwaveprj (iatm,2)%cp(1:2,jlmn) 830 cpj1(1:2,1)=cwaveprj1(iatm,1)%cp(1:2,jlmn) 831 cpj1(1:2,2)=cwaveprj1(iatm,2)%cp(1:2,jlmn) 832 do ilmn=1,jlmn 833 klmn=j0lmn+ilmn 834 klmn_re=cplex_rhoij*(klmn-1)+1 835 cpi0(1:2,1)=cwaveprj (iatm,1)%cp(1:2,ilmn) 836 cpi0(1:2,2)=cwaveprj (iatm,2)%cp(1:2,ilmn) 837 cpi1(1:2,1)=cwaveprj1(iatm,1)%cp(1:2,ilmn) 838 cpi1(1:2,2)=cwaveprj1(iatm,2)%cp(1:2,ilmn) 839 ro11_re=zero;ro22_re=zero 840 ro12_re=zero;ro21_re=zero 841 ro12_im=zero;ro21_im=zero 842 do iplex=1,cplex 843 ro11_re=ro11_re+cpj0(iplex,1)*cpi1(iplex,1)+cpi0(iplex,1)*cpj1(iplex,1) 844 ro22_re=ro22_re+cpj0(iplex,2)*cpi1(iplex,2)+cpi0(iplex,2)*cpj1(iplex,2) 845 end do 846 pawrhoij(iatom)%rhoij_(klmn_re,1)=pawrhoij(iatom)%rhoij_(klmn_re,1)+weight*(ro11_re+ro22_re) 847 if (nspden_rhoij>1) then 848 do iplex=1,cplex 849 ro12_re=ro12_re+cpj0(iplex,1)*cpi1(iplex,2)+cpi0(iplex,2)*cpj1(iplex,1) 850 ro21_re=ro21_re+cpj0(iplex,2)*cpi1(iplex,1)+cpi0(iplex,1)*cpj1(iplex,2) 851 end do 852 if (.not.use_zeromag) then 853 pawrhoij(iatom)%rhoij_(klmn_re,4)=pawrhoij(iatom)%rhoij_(klmn_re,4)+weight*(ro11_re-ro22_re) 854 pawrhoij(iatom)%rhoij_(klmn_re,2)=pawrhoij(iatom)%rhoij_(klmn_re,2)+weight*(ro12_re+ro21_re) 855 end if 856 if (cplex==2) then 857 !Important note: the present implementation follows eq(15) in Hobbs et al, PRB 62, 11556(2000) 858 ! rho^alpha,beta_ij = Sum[<Psi^beta|pi><pj|Psi^alpha] (alpha and beta exponents inverted) 859 ro12_im=cpj0(2,1)*cpi1(1,2)-cpi1(2,2)*cpj0(1,1)+cpi0(1,2)*cpj1(2,1)-cpj1(1,1)*cpi0(2,2) 860 ro21_im=cpj0(2,2)*cpi1(1,1)-cpi1(2,1)*cpj0(1,2)+cpi0(1,1)*cpj1(2,2)-cpj1(1,2)*cpi0(2,1) 861 if (.not.use_zeromag) then 862 pawrhoij(iatom)%rhoij_(klmn_re,3)=pawrhoij(iatom)%rhoij_(klmn_re,3)+weight*(ro21_im-ro12_im) 863 end if 864 end if 865 end if 866 if (compute_impart) then 867 klmn_im=klmn_re+iq0 868 if (nspden_rhoij>1) pawrhoij(iatom)%rhoij_(klmn_re,3)=pawrhoij(iatom)%rhoij_(klmn_re,3)+weight*(ro12_re-ro21_re) 869 if (cplex==2) then 870 ro11_im=cpj0(2,1)*cpi1(1,1)-cpi1(2,1)*cpj0(1,1)+cpi0(1,1)*cpj1(2,1)-cpj1(1,1)*cpi0(2,1) 871 ro22_im=cpj0(2,2)*cpi1(1,2)-cpi1(2,2)*cpj0(1,2)+cpi0(1,2)*cpj1(2,2)-cpj1(1,2)*cpi0(2,2) 872 pawrhoij(iatom)%rhoij_(klmn_im,1)=pawrhoij(iatom)%rhoij_(klmn_im,1)+weight*(ro11_im+ro22_im) 873 if (nspden_rhoij>1) then 874 pawrhoij(iatom)%rhoij_(klmn_im,4)=pawrhoij(iatom)%rhoij_(klmn_im,4)+weight*(ro11_im-ro22_im) 875 pawrhoij(iatom)%rhoij_(klmn_im,2)=pawrhoij(iatom)%rhoij_(klmn_re,2)+weight*(ro12_im+ro21_im) 876 end if 877 end if 878 end if 879 end do 880 end do 881 end do 882 end if 883 884 ! Accumulate (n,k) contribution to rhoij1 885 ! due to derivative of projectors 886 if (ipert/=natom+2) then 887 compute_impart=(pawrhoij(1)%cplex_rhoij==2) 888 compute_impart_cplex=(compute_impart.and.(cplex==2)) 889 substract_diagonal=(ipert==natom+3) 890 if (nspinor==1) then 891 do iatom=1,my_natom 892 iatom1=iatom;if (paral_atom) iatom1=my_atmtab(iatom) 893 iatm=atindx(iatom1) 894 if (ipert<=natom.and.iatom/=ipert) cycle 895 cplex_rhoij=pawrhoij(iatom)%cplex_rhoij 896 do jlmn=1,pawrhoij(iatom)%lmn_size 897 j0lmn=jlmn*(jlmn-1)/2 898 cpj0 (1:2,1) =cwaveprj(iatm,1)%cp (1:2 ,jlmn) 899 dcpj0(1:2,1,1)=cwaveprj(iatm,1)%dcp(1:2,1,jlmn) 900 do ilmn=1,jlmn 901 klmn=j0lmn+ilmn 902 klmn_re=cplex_rhoij*(klmn-1)+1 903 cpi0 (1:2,1) =cwaveprj(iatm,1)%cp (1:2 ,ilmn) 904 dcpi0(1:2,1,1)=cwaveprj(iatm,1)%dcp(1:2,1,ilmn) 905 ro11_re=zero 906 do iplex=1,cplex 907 ro11_re=ro11_re+dcpi0(iplex,1,1)*cpj0(iplex,1)+cpi0(iplex,1)*dcpj0(iplex,1,1) 908 end do 909 if (substract_diagonal) then 910 do iplex=1,cplex 911 ro11_re=ro11_re-cpi0(iplex,1)*cpj0(iplex,1) 912 end do 913 end if 914 pawrhoij(iatom)%rhoij_(klmn_re,isppol)=pawrhoij(iatom)%rhoij_(klmn_re,isppol)+weight*ro11_re 915 ! This imaginary part does not have to be computed 916 ! It is cancelled because rho_ij+rho_ji is stored in rho_ij 917 if (compute_impart_cplex) then 918 klmn_im=klmn_re+1 919 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) 920 if (substract_diagonal) then 921 ro11_im=ro11_im-cpi0(1,1)*cpj0(2,1)+cpi0(2,1)*cpj0(1,1) 922 end if 923 pawrhoij(iatom)%rhoij_(klmn_im,isppol)=pawrhoij(iatom)%rhoij_(klmn_im,isppol)+weight*ro11_im 924 end if 925 end do 926 end do 927 end do 928 else ! nspinor=2 929 do iatom=1,my_natom 930 iatom1=iatom;if (paral_atom) iatom1=my_atmtab(iatom) 931 iatm=atindx(iatom1) 932 if (ipert<=natom.and.iatom/=ipert) cycle 933 cplex_rhoij=pawrhoij(iatom)%cplex_rhoij 934 nspden_rhoij=pawrhoij(iatom)%nspden 935 do jlmn=1,pawrhoij(iatom)%lmn_size 936 j0lmn=jlmn*(jlmn-1)/2 937 cpj0 (1:2,1) =cwaveprj(iatm,1)%cp (1:2 ,jlmn) 938 dcpj0(1:2,1,1)=cwaveprj(iatm,1)%dcp(1:2,1,jlmn) 939 cpj0 (1:2,2) =cwaveprj(iatm,2)%cp (1:2 ,jlmn) 940 dcpj0(1:2,2,1)=cwaveprj(iatm,2)%dcp(1:2,1,jlmn) 941 do ilmn=1,jlmn 942 klmn=j0lmn+ilmn 943 klmn_re=cplex_rhoij*(klmn-1)+1 944 cpi0 (1:2,1) =cwaveprj(iatm,1)%cp (1:2 ,ilmn) 945 dcpi0(1:2,1,1)=cwaveprj(iatm,1)%dcp(1:2,1,ilmn) 946 cpi0 (1:2,2) =cwaveprj(iatm,2)%cp (1:2 ,ilmn) 947 dcpi0(1:2,2,1)=cwaveprj(iatm,2)%dcp(1:2,1,ilmn) 948 ro11_re=zero;ro22_re=zero 949 ro12_re=zero;ro21_re=zero 950 ro12_im=zero;ro21_im=zero 951 do iplex=1,cplex 952 ro11_re=ro11_re+dcpi0(iplex,1,1)*cpj0(iplex,1)+cpi0(iplex,1)*dcpj0(iplex,1,1) 953 ro22_re=ro22_re+dcpi0(iplex,2,1)*cpj0(iplex,2)+cpi0(iplex,2)*dcpj0(iplex,2,1) 954 end do 955 if (substract_diagonal) then 956 do iplex=1,cplex 957 ro11_re=ro11_re-cpi0(iplex,1)*cpj0(iplex,1) 958 ro22_re=ro22_re-cpi0(iplex,2)*cpj0(iplex,2) 959 end do 960 end if 961 pawrhoij(iatom)%rhoij_(klmn_re,1)=pawrhoij(iatom)%rhoij_(klmn_re,1)+weight*(ro11_re+ro22_re) 962 if (nspden_rhoij>1) then 963 do iplex=1,cplex 964 ro12_re=ro12_re+dcpi0(iplex,2,1)*cpj0(iplex,1)+cpi0(iplex,2)*dcpj0(iplex,1,1) 965 ro21_re=ro21_re+dcpi0(iplex,1,1)*cpj0(iplex,2)+cpi0(iplex,1)*dcpj0(iplex,2,1) 966 end do 967 if (substract_diagonal) then 968 do iplex=1,cplex 969 ro12_re=ro12_re-cpi0(iplex,2)*cpj0(iplex,1) 970 ro21_re=ro21_re-cpi0(iplex,1)*cpj0(iplex,2) 971 end do 972 end if 973 if (.not.use_zeromag) then 974 pawrhoij(iatom)%rhoij_(klmn_re,4)=pawrhoij(iatom)%rhoij_(klmn_re,4)+weight*(ro11_re-ro22_re) 975 pawrhoij(iatom)%rhoij_(klmn_re,2)=pawrhoij(iatom)%rhoij_(klmn_re,2)+weight*(ro12_re+ro21_re) 976 end if 977 if (cplex==2) then 978 !Important note: the present implementation follows eq(15) in Hobbs et al, PRB 62, 11556(2000) 979 ! rho^alpha,beta_ij = Sum[<Psi^beta|pi><pj|Psi^alpha] (alpha and beta exponents inverted) 980 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) 981 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) 982 if (substract_diagonal) then 983 ro12_im=ro12_im-cpi0(1,2)*cpj0(2,1)+cpi0(2,2)*cpj0(1,1) 984 ro21_im=ro21_im-cpi0(1,1)*cpj0(2,2)+cpi0(2,1)*cpj0(1,2) 985 end if 986 if (.not.use_zeromag) then 987 pawrhoij(iatom)%rhoij_(klmn_re,3)=pawrhoij(iatom)%rhoij_(klmn_re,3)+weight*(ro21_im-ro12_im) 988 end if 989 end if 990 end if 991 if (compute_impart) then 992 klmn_im=klmn_re+1 993 if (nspden_rhoij>1) pawrhoij(iatom)%rhoij_(klmn_im,3)=pawrhoij(iatom)%rhoij_(klmn_im,3)+weight*(ro12_re-ro21_re) 994 if (cplex==2) then 995 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) 996 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) 997 if (substract_diagonal) then 998 ro11_im=ro11_im-cpi0(1,1)*cpj0(2,1)+cpi0(2,1)*cpj0(1,1) 999 ro22_im=ro22_im-cpi0(1,2)*cpj0(2,2)+cpi0(2,2)*cpj0(1,2) 1000 end if 1001 pawrhoij(iatom)%rhoij_(klmn_im,1)=pawrhoij(iatom)%rhoij_(klmn_im,1)+weight*(ro11_im+ro22_im) 1002 if (nspden_rhoij>1) then 1003 pawrhoij(iatom)%rhoij_(klmn_im,4)=pawrhoij(iatom)%rhoij_(klmn_im,4)+weight*(ro11_im-ro22_im) 1004 pawrhoij(iatom)%rhoij_(klmn_im,2)=pawrhoij(iatom)%rhoij_(klmn_im,2)+weight*(ro12_im+ro21_im) 1005 end if 1006 end if 1007 end if 1008 end do 1009 end do 1010 end do 1011 end if 1012 end if 1013 1014 else if (option==3) then 1015 1016 ! ================================================================== 1017 ! === OPTION 3: Accumulate (n,k) contribution to drhoij/dr ========= 1018 ! ================================================================== 1019 1020 compute_impart=((.not.use_timerev).and.(pawrhoij(1)%cplex_rhoij==2)) 1021 compute_impart_cplex=((compute_impart).and.(cplex==2)) 1022 if (nspinor==1) then 1023 do iatom=1,my_natom 1024 iatom1=iatom;if (paral_atom) iatom1=my_atmtab(iatom) 1025 iatm=atindx(iatom1) 1026 cplex_rhoij=pawrhoij(iatom)%cplex_rhoij 1027 do jlmn=1,pawrhoij(iatom)%lmn_size 1028 j0lmn=jlmn*(jlmn-1)/2 1029 cpj0(1:cplex,1) =cwaveprj(iatm,1)%cp (1:cplex,jlmn) 1030 dcpj0(1:cplex,1,1:ncpgr)=cwaveprj(iatm,1)%dcp(1:cplex,1:ncpgr,jlmn) 1031 do ilmn=1,jlmn 1032 klmn=j0lmn+ilmn 1033 klmn_re=cplex_rhoij*(klmn-1)+1 1034 cpi0(1:cplex,1) =cwaveprj(iatm,1)%cp (1:cplex,ilmn) 1035 dcpi0(1:cplex,1,1:ncpgr)=cwaveprj(iatm,1)%dcp(1:cplex,1:ncpgr,ilmn) 1036 do mu=1,ncpgr 1037 ro11_re=zero 1038 do iplex=1,cplex 1039 ro11_re=ro11_re+dcpi0(iplex,1,mu)*cpj0(iplex,1)+cpi0(iplex,1)*dcpj0(iplex,1,mu) 1040 end do 1041 pawrhoij(iatom)%grhoij(mu,klmn_re,isppol)=pawrhoij(iatom)%grhoij(mu,klmn_re,isppol)+weight*ro11_re 1042 end do 1043 if (compute_impart_cplex) then 1044 klmn_im=klmn_re+1 1045 do mu=1,ncpgr 1046 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) 1047 pawrhoij(iatom)%grhoij(mu,klmn_im,isppol)=pawrhoij(iatom)%grhoij(mu,klmn_im,isppol)+weight*ro11_im 1048 end do 1049 end if 1050 end do 1051 end do 1052 end do 1053 else ! nspinor=2 1054 do iatom=1,my_natom 1055 iatom1=iatom;if (paral_atom) iatom1=my_atmtab(iatom) 1056 iatm=atindx(iatom1) 1057 cplex_rhoij=pawrhoij(iatom)%cplex_rhoij 1058 nspden_rhoij=pawrhoij(iatom)%nspden 1059 do jlmn=1,pawrhoij(iatom)%lmn_size 1060 j0lmn=jlmn*(jlmn-1)/2 1061 cpj0(1:cplex,1) =cwaveprj(iatm,1)%cp (1:cplex,jlmn) 1062 cpj0(1:cplex,2) =cwaveprj(iatm,2)%cp (1:cplex,jlmn) 1063 dcpj0(1:cplex,1,1:ncpgr)=cwaveprj(iatm,1)%dcp(1:cplex,1:ncpgr,jlmn) 1064 dcpj0(1:cplex,2,1:ncpgr)=cwaveprj(iatm,2)%dcp(1:cplex,1:ncpgr,jlmn) 1065 do ilmn=1,jlmn 1066 klmn=j0lmn+ilmn 1067 klmn_re=cplex_rhoij*(klmn-1)+1 1068 klmn_im=klmn_re+1 1069 cpi0(1:cplex,1) =cwaveprj(iatm,1)%cp (1:cplex,ilmn) 1070 cpi0(1:cplex,2) =cwaveprj(iatm,2)%cp (1:cplex,ilmn) 1071 dcpi0(1:cplex,1,1:ncpgr)=cwaveprj(iatm,1)%dcp(1:cplex,1:ncpgr,ilmn) 1072 dcpi0(1:cplex,2,1:ncpgr)=cwaveprj(iatm,2)%dcp(1:cplex,1:ncpgr,ilmn) 1073 do mu=1,ncpgr 1074 ro11_re=zero;ro22_re=zero 1075 ro12_re=zero;ro21_re=zero 1076 ro12_im=zero;ro21_im=zero 1077 do iplex=1,cplex 1078 ro11_re=ro11_re+dcpi0(iplex,1,mu)*cpj0(iplex,1)+cpi0(iplex,1)*dcpj0(iplex,1,mu) 1079 ro22_re=ro22_re+dcpi0(iplex,2,mu)*cpj0(iplex,2)+cpi0(iplex,2)*dcpj0(iplex,2,mu) 1080 end do 1081 pawrhoij(iatom)%grhoij(mu,klmn_re,1)=pawrhoij(iatom)%grhoij(mu,klmn_re,1)+weight*(ro11_re+ro22_re) 1082 if (nspden_rhoij>1) then 1083 do iplex=1,cplex 1084 ro12_re=ro12_re+dcpi0(iplex,2,mu)*cpj0(iplex,1)+cpi0(iplex,2)*dcpj0(iplex,1,mu) 1085 ro21_re=ro21_re+dcpi0(iplex,1,mu)*cpj0(iplex,2)+cpi0(iplex,1)*dcpj0(iplex,2,mu) 1086 end do 1087 if (.not.use_zeromag) then 1088 pawrhoij(iatom)%grhoij(mu,klmn_re,4)=pawrhoij(iatom)%grhoij(mu,klmn_re,4)+weight*(ro11_re-ro22_re) 1089 pawrhoij(iatom)%grhoij(mu,klmn_re,2)=pawrhoij(iatom)%grhoij(mu,klmn_re,2)+weight*(ro12_re+ro21_re) 1090 end if 1091 if (cplex==2) then 1092 !Important note: the present implementation follows eq(15) in Hobbs et al, PRB 62, 11556(2000) 1093 ! rho^alpha,beta_ij = Sum[<Psi^beta|pi><pj|Psi^alpha] (alpha and beta exponents inverted) 1094 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) 1095 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) 1096 if (.not.use_zeromag) then 1097 pawrhoij(iatom)%grhoij(mu,klmn_re,3)=pawrhoij(iatom)%grhoij(mu,klmn_re,3)+weight*(ro21_im-ro12_im) 1098 end if 1099 end if 1100 end if 1101 if (compute_impart) then 1102 if (nspden_rhoij>1) then 1103 pawrhoij(iatom)%grhoij(mu,klmn_im,3)=pawrhoij(iatom)%grhoij(mu,klmn_im,3)+weight*(ro12_re-ro21_re) 1104 end if 1105 if (cplex==2) then 1106 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) 1107 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) 1108 pawrhoij(iatom)%grhoij(mu,klmn_im,1)=pawrhoij(iatom)%grhoij(mu,klmn_im,1)+weight*(ro11_im+ro22_im) 1109 if (nspden_rhoij>1) then 1110 pawrhoij(iatom)%grhoij(mu,klmn_im,4)=pawrhoij(iatom)%grhoij(mu,klmn_im,4)+weight*(ro11_im-ro22_im) 1111 pawrhoij(iatom)%grhoij(mu,klmn_im,2)=pawrhoij(iatom)%grhoij(mu,klmn_im,2)+weight*(ro12_im+ro21_im) 1112 end if 1113 end if 1114 end if 1115 end do 1116 end do 1117 end do 1118 end do 1119 end if 1120 1121 ! End 1122 end if ! option 1123 1124 !Destroy atom table used for parallelism 1125 call free_my_atmtab(my_atmtab,my_atmtab_allocated) 1126 1127 DBG_EXIT("COLL") 1128 1129 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
nspden=number of spin components for the density
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
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)
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
102 subroutine pawmkrhoij(atindx,atindx1,cprj,dimcprj,istwfk,kptopt,mband,mband_cprj,mcprj,mkmem,mpi_enreg,& 103 & natom,nband,nkpt,nspden,nspinor,nsppol,occ,paral_kgb,paw_dmft,pawrhoij,unpaw,usewvl,wtk) 104 !Arguments --------------------------------------------- 105 !scalars 106 integer,intent(in) :: kptopt,mband,mband_cprj,mcprj,mkmem,natom,nkpt,nspden,nspinor,nsppol 107 integer,intent(in) :: paral_kgb,unpaw,usewvl 108 type(MPI_type),intent(in) :: mpi_enreg 109 !arrays 110 integer,intent(in) :: atindx(natom),atindx1(natom),dimcprj(natom),istwfk(nkpt) 111 integer,intent(in) :: nband(nkpt*nsppol) 112 real(dp),intent(in) :: occ(mband*nkpt*nsppol),wtk(nkpt) 113 type(pawcprj_type),target,intent(in) :: cprj(natom,mcprj) 114 type(paw_dmft_type),intent(in) :: paw_dmft 115 type(pawrhoij_type),intent(inout),target:: pawrhoij(:) 116 117 !Local variables --------------------------------------- 118 !scalars 119 integer,parameter :: max_nband_cprj=100 120 integer :: bdtot_index,cplex 121 integer :: cplex_rhoij,iatom,iatom_tot,ib,ib1,iband,ibc1,ibg,ib_this_proc,ierr 122 integer :: ikpt,iorder_cprj,isppol,jb_this_proc,jbg,me,my_nspinor,nband_k,nband_k_cprj 123 integer :: nbandc1,nband_k_cprj_read,nband_k_cprj_used,nprocband,nrhoij,nspden_rhoij 124 integer :: option,spaceComm,use_nondiag_occup_dmft 125 logical :: locc_test,paral_atom,use_timerev,use_zeromag 126 integer :: ib1_this_proc, ib_loop, proc_sender, proc_recver 127 real(dp) :: wtk_k 128 character(len=500) :: msg 129 130 !arrays 131 integer :: n2buff 132 integer, allocatable :: typat(:),req_correl(:,:,:) 133 real(dp) :: occup(2) 134 real(dp) ABI_ASYNC, allocatable :: buffer_cprj_correl(:,:,:) 135 character(len=8),parameter :: dspin(6)=(/"up ","down ","dens (n)","magn (x)","magn (y)","magn (z)"/) 136 type(pawcprj_type),allocatable :: cprj_tmp(:,:),cwaveprj(:,:),cwaveprjb(:,:) 137 type(pawcprj_type),pointer :: cprj_ptr(:,:) 138 type(pawrhoij_type),pointer :: pawrhoij_all(:) 139 140 !************************************************************************ 141 142 DBG_ENTER("COLL") 143 144 ABI_CHECK(mkmem/=0,"mkmem==0 not supported anymore!") 145 146 !Init MPI data 147 ! spaceComm=mpi_enreg%comm_cell 148 ! if(mpi_enreg%paral_kgb==1) spaceComm=mpi_enreg%comm_kpt 149 spaceComm=mpi_enreg%comm_kpt 150 me=mpi_enreg%me_kpt 151 152 !Check size of cprj 153 my_nspinor=max(1,nspinor/mpi_enreg%nproc_spinor) 154 if (mcprj/=my_nspinor*mband_cprj*mkmem*nsppol) then 155 msg=' wrong size for cprj !' 156 ABI_BUG(msg) 157 end if 158 159 !Check if cprj is distributed over bands 160 nprocband=(mband/mband_cprj) 161 if (paral_kgb==1.and.nprocband/=mpi_enreg%nproc_band) then 162 msg='mband/mband_cprj must be equal to nproc_band!' 163 ABI_BUG(msg) 164 end if 165 166 if( usewvl==1 .and. (nprocband/=1)) then 167 write(msg,'(2a)') ch10,& 168 & 'Parallelization over bands is not compatible with WAVELETS!' 169 ABI_ERROR(msg) 170 end if 171 172 !Initialise and check dmft variables 173 if(paw_dmft%use_sc_dmft/=0) then 174 nbandc1=paw_dmft%mbandc 175 else 176 nbandc1=1 177 end if 178 179 !Size of pawrhoij datastructure 180 nrhoij=size(pawrhoij) 181 182 !Check if pawrhoij is distributed over atomic sites 183 paral_atom=(nrhoij/=natom.and.mpi_enreg%nproc_atom>1) 184 if (paral_atom.and.nrhoij/=mpi_enreg%my_natom) then 185 msg='Size of pawrhoij should be natom or my_natom!' 186 ABI_BUG(msg) 187 end if 188 189 !Allocate temporary cwaveprj storage 190 ABI_MALLOC(cwaveprj,(natom,nspinor)) 191 call pawcprj_alloc(cwaveprj,0,dimcprj) 192 if(paw_dmft%use_sc_dmft/=0) then 193 ABI_MALLOC(cwaveprjb,(natom,nspinor)) 194 call pawcprj_alloc(cwaveprjb,0,dimcprj) 195 end if 196 197 if (paw_dmft%use_sc_dmft /= 0 .and. mpi_enreg%paral_kgb /= 0) then 198 if(paw_dmft%use_bandc(mpi_enreg%me_band+1)) then 199 n2buff = nspinor*sum(dimcprj) 200 ABI_MALLOC(buffer_cprj_correl,(2,n2buff,nbandc1)) 201 ABI_MALLOC(req_correl,(nbandc1, nkpt, nsppol)) 202 req_correl(:,:,:) = 0 203 end if 204 end if 205 206 !Initialize temporary file (if used) 207 iorder_cprj=0 208 209 !Build and initialize unpacked rhoij (to be computed here) 210 call pawrhoij_init_unpacked(pawrhoij) 211 212 !If pawrhoij is MPI-distributed over atomic sites, we have to gather it 213 if (paral_atom) then 214 ABI_MALLOC(pawrhoij_all,(natom)) 215 ABI_MALLOC(typat,(natom)) 216 typat(:)=0;cplex_rhoij=0;nspden_rhoij=0 217 do iatom=1,nrhoij 218 iatom_tot=mpi_enreg%my_atmtab(iatom) 219 typat(iatom_tot)=pawrhoij(iatom)%itypat 220 cplex_rhoij=max(cplex_rhoij,pawrhoij(iatom)%cplex_rhoij) 221 nspden_rhoij=max(nspden_rhoij,pawrhoij(iatom)%nspden) 222 end do 223 call xmpi_sum(typat,mpi_enreg%comm_atom,ierr) 224 call xmpi_max_ip(cplex_rhoij,mpi_enreg%comm_atom,ierr) 225 call xmpi_max_ip(nspden_rhoij,mpi_enreg%comm_atom,ierr) 226 call pawrhoij_alloc(pawrhoij_all,cplex_rhoij,nspden_rhoij,nspinor,& 227 & nsppol,typat,lmnsize=dimcprj,use_rhoijp=0,use_rhoij_=1) 228 ABI_FREE(typat) 229 else 230 pawrhoij_all => pawrhoij 231 end if 232 233 !LOOP OVER SPINS 234 option=1 235 use_timerev=(kptopt>0.and.kptopt<3) 236 use_zeromag=(pawrhoij_all(1)%nspden==4.and.nspden==1) 237 bdtot_index=0;ibg=0;jbg=0 238 do isppol=1,nsppol 239 240 ! LOOP OVER k POINTS 241 do ikpt=1,nkpt 242 243 nband_k=nband(ikpt+(isppol-1)*nkpt) 244 nband_k_cprj=nband_k/nprocband 245 wtk_k=wtk(ikpt) 246 247 if(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,isppol,me)) then 248 bdtot_index=bdtot_index+nband_k 249 cycle 250 end if 251 252 cplex=2;if (istwfk(ikpt)>1) cplex=1 253 254 ! In case of spinors parallelism, need some extra storage 255 if (mpi_enreg%paral_spinor==1) then 256 nband_k_cprj_used=min(max_nband_cprj,nband_k_cprj) 257 ABI_MALLOC(cprj_tmp,(natom,my_nspinor*nband_k_cprj_used)) 258 ABI_MALLOC(cprj_ptr,(natom, nspinor*nband_k_cprj_used)) 259 call pawcprj_alloc(cprj_tmp,0,dimcprj) 260 call pawcprj_alloc(cprj_ptr,0,dimcprj) 261 else 262 cprj_ptr => cprj 263 end if 264 265 ! In case of band parallelism combined with self consistent DMFT, need to 266 ! exchange bands cprj 267 if (paw_dmft%use_sc_dmft /= 0 .and. mpi_enreg%paral_kgb /= 0) then 268 if (paw_dmft%use_bandc(mpi_enreg%me_band+1)) then 269 ! only proc using correlated band have to do this 270 do ibc1=1,nbandc1 271 proc_sender = paw_dmft%bandc_proc(ibc1) 272 if(proc_sender == mpi_enreg%me_band) then 273 274 ! get the index of band local to this proc 275 ib1 = paw_dmft%include_bands(ibc1) 276 ib1_this_proc = 0 277 do ib_loop=1,ib1-1 278 if (mod((ib_loop-1)/mpi_enreg%bandpp,mpi_enreg%nproc_band) == mpi_enreg%me_band) then 279 ib1_this_proc = ib1_this_proc+1 280 end if 281 end do 282 ib1_this_proc = ib1_this_proc+1 283 284 ! extract the band 285 call pawcprj_get(atindx1,cwaveprjb,cprj_ptr,natom,ib1_this_proc,ibg,ikpt,& 286 & iorder_cprj,isppol,mband_cprj,mkmem,natom,1,nband_k_cprj,nspinor,nsppol,& 287 & unpaw,mpicomm=mpi_enreg%comm_kpt,proc_distrb=mpi_enreg%proc_distrb) 288 289 call pawcprj_pack(dimcprj,cwaveprjb,buffer_cprj_correl(:,:,ibc1)) 290 do proc_recver=0,mpi_enreg%nproc_band-1 291 if (proc_sender /= proc_recver .and. paw_dmft%use_bandc(proc_recver+1)) then 292 ! locc_test = At least one of the bands used by proc_recver have a non neglectable occnd 293 locc_test = .false. 294 do ib_loop=1,nbandc1 295 if(proc_recver == paw_dmft%bandc_proc(ib_loop)) then 296 ib = paw_dmft%include_bands(ib_loop) 297 locc_test = locc_test .or. (abs(paw_dmft%occnd(1,ib,ib1,ikpt,isppol))+& 298 & abs(paw_dmft%occnd(2,ib,ib1,ikpt,isppol))>tol8) 299 end if 300 end do 301 if(locc_test) then 302 ! send to proc_recver 303 ierr = 0 304 call xmpi_isend(buffer_cprj_correl(:,:,ibc1),proc_recver,& 305 & 10000+ibc1+nbandc1*(ikpt+nsppol*isppol),mpi_enreg%comm_band,& 306 & req_correl(ibc1,ikpt,isppol),ierr) 307 ! force sending or buffering 308 call xmpi_wait(req_correl(ibc1,ikpt,isppol), ierr) 309 end if 310 end if 311 end do 312 else 313 ! locc_test = At least one of the bands used by this proc have a non neglectable occnd 314 locc_test = .false. 315 do ib_loop=1,nbandc1 316 if(mpi_enreg%me_band == paw_dmft%bandc_proc(ib_loop)) then 317 ib = paw_dmft%include_bands(ib_loop) 318 ib1 = paw_dmft%include_bands(ibc1) 319 locc_test = locc_test .or. (abs(paw_dmft%occnd(1,ib,ib1,ikpt,isppol))+& 320 & abs(paw_dmft%occnd(2,ib,ib1,ikpt,isppol))>tol8) 321 end if 322 end do 323 if(locc_test) then 324 ! recv from proc_sender 325 ierr = 0 326 call xmpi_irecv(buffer_cprj_correl(:,:,ibc1),proc_sender,& 327 & 10000+ibc1+nbandc1*(ikpt+nsppol*isppol),mpi_enreg%comm_band,& 328 & req_correl(ibc1,ikpt,isppol),ierr) 329 end if 330 end if 331 end do 332 end if 333 end if 334 335 ierr = 0 336 ! LOOP OVER BANDS 337 ib_this_proc=0;jb_this_proc=0 338 do ib=1,nband_k 339 iband=bdtot_index+ib 340 341 ! Parallelization: treat only some bands 342 if(xmpi_paral==1)then 343 if (paral_kgb==1) then 344 if (mod((ib-1)/mpi_enreg%bandpp,mpi_enreg%nproc_band)/=mpi_enreg%me_band) cycle 345 else 346 if (mpi_enreg%proc_distrb(ikpt,ib,isppol)/=me) cycle 347 end if 348 end if 349 ib_this_proc=ib_this_proc+1 350 351 ! In case of spinors parallelism, gather cprj because we need both components together 352 ! We do that nband_k_cprj_used by nband_k_cprj_used bands 353 if (mpi_enreg%paral_spinor==1) then 354 jb_this_proc=jb_this_proc+1 355 if (mod(jb_this_proc,nband_k_cprj_used)==1) then 356 ib_this_proc=1 357 nband_k_cprj_read=nband_k_cprj_used 358 if (nband_k_cprj<jb_this_proc+nband_k_cprj_used-1) nband_k_cprj_read=nband_k_cprj-jb_this_proc+1 359 call pawcprj_get(atindx1,cprj_tmp,cprj,natom,jb_this_proc,jbg,ikpt,iorder_cprj,isppol,& 360 & mband_cprj,mkmem,natom,nband_k_cprj_read,nband_k_cprj,my_nspinor,nsppol,unpaw,& 361 & mpicomm=mpi_enreg%comm_kpt,proc_distrb=mpi_enreg%proc_distrb) 362 call pawcprj_gather_spin(cprj_tmp,cprj_ptr,natom,nband_k_cprj_read,my_nspinor,nspinor,& 363 & mpi_enreg%comm_spinor,ierr) 364 end if 365 end if 366 367 ! DMFT: LOOP ON ADDITIONAL BANDS 368 do ibc1=1,nbandc1 369 ! check if dmft and occupations 370 ! write(std_out,*) 'ib,ibc1 ',ib,ibc1 371 372 ! if ib is not part a band correlated in dmft do not repeat the following 373 if(ibc1 /= 1) then 374 if (paw_dmft%use_sc_dmft == 0) then 375 cycle 376 else 377 if (.not.(paw_dmft%band_in(ib))) cycle 378 end if 379 end if 380 381 ! DMFT stuff: extract cprj and occupations for additional band 382 if(paw_dmft%use_sc_dmft /= 0) then 383 if(paw_dmft%band_in(ib)) then 384 ! write(std_out,*) 'use_sc_dmft=1 ib,ib1',ib,ib1 385 ! write(std_out,*) 'ib, ib1 ',paw_dmft%band_in(ib),paw_dmft%band_in(ib1) 386 387 ib1 = paw_dmft%include_bands(ibc1) ! indice reel de la bande 388 389 use_nondiag_occup_dmft = 1 390 locc_test = abs(paw_dmft%occnd(1,ib,ib1,ikpt,isppol))+abs(paw_dmft%occnd(2,ib,ib1,ikpt,isppol))>tol8 391 392 occup(1) = paw_dmft%occnd(1,ib,ib1,ikpt,isppol) 393 occup(2) = paw_dmft%occnd(2,ib,ib1,ikpt,isppol) 394 395 ! write(std_out,*) 'use_sc_dmft=1,band_in(ib)=1, ib,ibc1',ib,ib1,locc_test 396 ! 397 if (locc_test .or. mkmem == 0) then 398 399 if (paral_kgb==1) then 400 ! cprj have already been extracted 401 if (paw_dmft%bandc_proc(ibc1) /= mpi_enreg%me_band) then 402 ! if the band is not on this proc, wait for the recv to complete 403 ierr = 0 404 call xmpi_wait(req_correl(ibc1,ikpt,isppol), ierr) 405 end if 406 call pawcprj_unpack(dimcprj,cwaveprjb,buffer_cprj_correl(:,:,ibc1)) 407 else ! paral_kgb /= 0 408 call pawcprj_get(atindx1,cwaveprjb,cprj_ptr,natom,ib1,ibg,ikpt,iorder_cprj,isppol,& 409 & mband_cprj,mkmem,natom,1,nband_k_cprj,nspinor,nsppol,unpaw,& 410 & mpicomm=mpi_enreg%comm_kpt,proc_distrb=mpi_enreg%proc_distrb) 411 end if 412 end if 413 else ! nbandc1=1 414 use_nondiag_occup_dmft=0 415 locc_test = (abs(occ(iband))>tol8) 416 occup(1) = occ(iband) 417 end if 418 else ! nbandc1=1 419 use_nondiag_occup_dmft=0 420 locc_test = (abs(occ(iband))>tol8) 421 occup(1) = occ(iband) 422 end if 423 424 ! Extract cprj for current band 425 ! Must read cprj when mkmem=0 (even if unused) to have right pointer inside _PAW file 426 if (locc_test.or.mkmem==0) then 427 call pawcprj_get(atindx1,cwaveprj,cprj_ptr,natom,ib_this_proc,ibg,ikpt,iorder_cprj,isppol,& 428 & mband_cprj,mkmem,natom,1,nband_k_cprj,nspinor,nsppol,unpaw,& 429 & mpicomm=mpi_enreg%comm_kpt,proc_distrb=mpi_enreg%proc_distrb) 430 end if 431 432 ! Accumulate contribution from (occupied) current band 433 if (locc_test) then 434 if(use_nondiag_occup_dmft == 1) then 435 call pawaccrhoij(atindx,cplex,cwaveprj,cwaveprjb,0,isppol,nrhoij,natom,& 436 & nspinor,occup(1),option,pawrhoij_all,use_timerev,use_zeromag,wtk_k,& 437 & occ_k_2=occup(2)) 438 else 439 call pawaccrhoij(atindx,cplex,cwaveprj,cwaveprj ,0,isppol,nrhoij,natom,& 440 & nspinor,occup(1),option,pawrhoij_all,use_timerev,use_zeromag,wtk_k) 441 end if 442 end if 443 end do ! ib1c 444 end do ! ib 445 446 if (mpi_enreg%paral_spinor==1) then 447 call pawcprj_free(cprj_tmp) 448 call pawcprj_free(cprj_ptr) 449 ABI_FREE(cprj_tmp) 450 ABI_FREE(cprj_ptr) 451 else 452 nullify(cprj_ptr) 453 end if 454 455 bdtot_index=bdtot_index+nband_k 456 if (mkmem/=0) then 457 if (mpi_enreg%paral_spinor==0) then 458 ibg=ibg+ nspinor*nband_k_cprj 459 else 460 jbg=jbg+my_nspinor*nband_k_cprj 461 end if 462 end if 463 464 end do ! ikpt 465 end do ! isppol 466 467 !call xmpi_barrier(mpi_enreg%comm_band) 468 469 !deallocate temporary cwaveprj/cprj storage 470 call pawcprj_free(cwaveprj) 471 ABI_FREE(cwaveprj) 472 473 if(paw_dmft%use_sc_dmft/=0) then 474 call pawcprj_free(cwaveprjb) 475 ABI_FREE(cwaveprjb) 476 end if 477 478 if (allocated(buffer_cprj_correl)) then 479 ABI_FREE(buffer_cprj_correl) 480 ABI_FREE(req_correl) 481 end if 482 483 !MPI: need to exchange rhoij_ between procs 484 if (paral_kgb==1.and.nprocband>1) then 485 call pawrhoij_mpisum_unpacked(pawrhoij_all,spaceComm,comm2=mpi_enreg%comm_band) 486 else 487 call pawrhoij_mpisum_unpacked(pawrhoij_all,spaceComm) 488 end if 489 490 !In case of distribution over atomic sites, dispatch rhoij 491 if (paral_atom) then 492 do iatom=1,nrhoij 493 iatom_tot=mpi_enreg%my_atmtab(iatom) 494 pawrhoij(iatom)%rhoij_(:,:)=pawrhoij_all(iatom_tot)%rhoij_(:,:) 495 end do 496 call pawrhoij_free(pawrhoij_all) 497 ABI_FREE(pawrhoij_all) 498 end if 499 500 DBG_EXIT("COLL") 501 502 end subroutine pawmkrhoij