TABLE OF CONTENTS
ABINIT/dfpt_vtorho [ Functions ]
NAME
dfpt_vtorho
FUNCTION
This routine compute the new 1-density from a fixed 1-potential (vtrial1) but might also simply compute eigenvectors and eigenvalues. The main part of it is a wf update over all k points
COPYRIGHT
Copyright (C) 1998-2018 ABINIT group (DCA, XG, GMR, AR, DRH, MB, XW, 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
cg(2,mpw*nspinor*mband*mkmem*nsppol)=planewave coefficients of wavefunctions cgq(2,mpw1*nspinor*mband*mkqmem*nsppol)=pw coefficients of GS wavefunctions at k+q. cg1(2,mpw1*nspinor*mband*mk1mem*nsppol)=pw coefficients of RF wavefunctions at k,q. cplex: if 1, real space 1-order functions on FFT grid are REAL; if 2, COMPLEX cprj(natom,nspinor*mband*mkmem*nsppol*usecprj)= wave functions at k projected with non-local projectors: cprj=<p_i|Cnk> cprjq(natom,nspinor*mband*mkqmem*nsppol*usecprj)= wave functions at k+q projected with non-local projectors: cprjq=<p_i|Cnk+q> dbl_nnsclo=if 1, will double the value of dtset%nnsclo doccde_rbz(mband*nkpt_rbz*nsppol)=derivative of occ_rbz wrt the energy docckqde(mband*nkpt_rbz*nsppol)=derivative of occkq wrt the energy dtefield = variables related to response Berry-phase calculation dtfil <type(datafiles_type)>=variables related to files dtset <type(dataset_type)>=all input variables for this dataset eigenq(mband*nkpt_rbz*nsppol)=GS eigenvalues at k+q (hartree) eigen0(mband*nkpt_rbz*nsppol)=GS eigenvalues at k (hartree) fermie1=derivative of fermi energy wrt (strain) perturbation gmet(3,3)=reciprocal space metric tensor in bohr**-2. gprimd(3,3)=dimensional reciprocal space primitive translations idir=direction of the perturbation indsy1(4,nsym1,natom)=indirect indexing array for atom labels ipert=type of the perturbation irrzon1(nfft**(1-1/nsym1),2,(nspden/nsppol)-3*(nspden/4))=irreducible zone data istwfk_rbz(nkpt_rbz)=input option parameter that describes the storage of wfs kg(3,mpw*mkmem)=reduced planewave coordinates. kg1(3,mpw1*mk1mem)=reduced planewave coordinates at k+q, with RF k points kpt_rbz(3,nkpt_rbz)=reduced coordinates of k points. mband=maximum number of bands mkmem =number of k points treated by this node (GS data). mkqmem =number of k+q points treated by this node (GS data) mk1mem =number of k points treated by this node (RF data) mpw=maximum dimensioned size of npw or wfs at k mpw1=maximum dimensioned size of npw for wfs at k+q (also for 1-order wfs). my_natom=number of atoms treated by current processor natom=number of atoms in cell. nband_rbz(nkpt_rbz*nsppol)=number of bands at each RF k point for each spin ncpgr=number of gradients stored in cprj array (cprj=<p_i|Cnk>) nfftf= -PAW ONLY- number of FFT grid points for the fine grid (nfftf=nfft for norm-conserving potential runs - see comment in respfn.F90) nkpt_rbz=number of k points in the IBZ for this perturbation mpi_enreg=information about MPI parallelization npwarr(nkpt_rbz)=number of planewaves in basis at this GS k point npwar1(nkpt_rbz)=number of planewaves in basis at this RF k+q point nspden=number of spin-density components nsppol=1 for unpolarized, 2 for spin-polarized nsym1=number of symmetry elements in space group consistent with perturbation ntypat=number of types of atoms in unit cell. occkq(mband*nkpt_rbz*nsppol)=occupation number for each band (often 2) at each k+q point of the reduced Brillouin zone. occ_rbz(mband*nkpt_rbz*nsppol)=occupation number for each band and k (usually 2) optres=0: the new value of the density is computed in place of the input value 1: only the density residual is computed ; the input density is kept paw_ij(natom*usepaw) <type(paw_ij_type)>=paw arrays given on (i,j) channels for the GS paw_ij1(natom) <type(paw_ij_type)>=1st-order paw arrays given on (i,j) channels pawang <type(pawang_type)>=paw angular mesh and related data pawang1 <type(pawang_type)>=pawang datastructure containing only the symmetries preserving the perturbation pawfgr <type(pawfgr_type)>=fine grid parameters and related data pawfgrtab(natom*usepaw) <type(pawfgrtab_type)>=atomic data given on fine rectangular grid for the GS pawrad(ntypat*usepaw) <type(pawrad_type)>=paw radial mesh and related data pawrhoij(natom) <type(pawrhoij_type)>= paw rhoij occupancies and related data for the GS pawtab(ntypat*usepaw) <type(pawtab_type)>=paw tabulated starting data phnons1(2,dtset%nfft**(1-1/nsym1),(nspden/nsppol)-3*(nspden/4))=nonsymmorphic translation phases ph1d(2,3*(2*dtset%mgfft+1)*natom)=one-dimensional structure factor information prtvol=control print volume and debugging output psps <type(pseudopotential_type)>=variables related to pseudopotentials pwindall(max(mpw,mpw1)*mkmem,8,3) = array used to compute the overlap matrices qmat(2,dtefield%mband_occ,dtefield%mband_occ,nkpt,2,3) = inverse of the overlap matrix rmet(3,3)=real space metric (bohr**2) rprimd(3,3)=dimensional real space primitive translations symaf1(nsym1)=(anti)ferromagnetic part of symmetry operations symrc1(3,3,nsym1)=symmetry operations in reciprocal space symrl1(3,3,nsym1)=symmetry operations in real space ucvol=unit cell volume in bohr**3. usecprj= 1 if cprj, cprjq, cprj1 arrays are stored in memory useylmgr1= 1 if ylmgr1 array is allocated ddk<wfk_t)=struct info DDK file vtrial(nfftf,nspden)=GS Vtrial(r). vtrial1(cplex*nfftf,nspden)=INPUT RF Vtrial(r). wtk_rbz(nkpt_rbz)=weight assigned to each k point. xred(3,natom)=reduced dimensionless atomic coordinates ylm(mpw*mkmem,mpsang*mpsang*useylm)= real spherical harmonics for each G and k point ylm1(mpw1*mk1mem,mpsang*mpsang*useylm)= real spherical harmonics for each G and k+g point ylmgr1(mpw*mkmem,3,mpsang*mpsang*useylm)= gradients of real spherical harmonics for each G and k+g point
OUTPUT
cg1(2,mpw*nspinor*mband*mk1mem*nsppol)=updated wavefunctions, orthogonalized to the occupied states cg1_active(2,mpw1*nspinor*mband*mk1mem*nsppol)=pw coefficients of RF wavefunctions at k,q. They are orthogonalized to the active. eigen1(2*mband*mband*nkpt_rbz*nsppol)=array for holding eigenvalues (hartree) edocc=correction to 2nd-order total energy coming from changes of occupation eeig0=0th-order eigenenergies part of 2nd-order total energy ek0=0th-order kinetic energy part of 2nd-order total energy. ek1=1st-order kinetic energy part of 2nd-order total energy (not for phonons) eloc0=0th-order local (psp+vxc+Hart) part of 2nd-order total energy enl0=0th-order nonlocal pseudopot. part of 2nd-order total energy. enl1=1st-order nonlocal pseudopot. part of 2nd-order total energy. gh1c_set(2,mpw1*nspinor*mband*mk1mem*nsppol*dim_eig2rf)= set of <G|H^{(1)}|nK> gh0c1_set(2,mpw1*nspinor*mband*mk1mem*nsppol*dim_eig2rf)= set of <G|H^{(0)}|\Psi^{(1)}> The wavefunction is orthogonal to the active space (for metals). It is not coherent with cg1. resid(mband*nkpt_rbz*nsppol)=residuals for each band over all k points. residm=maximum value from resid array (except for nbdbuf highest bands) rhog1(2,nfftf)=RF electron density in reciprocal space ==== if optres==1 nres2=square of the norm of the residual nvresid1(cplex*nfftf,nspden)=1st-order density residual ==== if psps%usepaw==1 cprj1(natom,nspinor*mband*mk1mem*nsppol*usecprj)= 1st-order wave functions at k,q projected with non-local projectors: cprj1=<p_i|C1nk,q> where p_i is a non-local projector nhat1(cplex*nfftf,nspden*psps%usepaw)=1st-order compensation charge density
SIDE EFFECTS
pawrhoij1(natom) <type(pawrhoij_type)>= 1st-order paw rhoij occupancies and related data rhor1(cplex*nfftf,nspden)=RF electron density in electrons/bohr**3.
PARENTS
dfpt_scfcv
CHILDREN
destroy_hamiltonian,destroy_rf_hamiltonian,dfpt_vtowfk,dfptff_gbefd dfptff_gradberry,fftpac,getgh1c_setup,init_hamiltonian init_rf_hamiltonian,load_spin_hamiltonian,load_spin_rf_hamiltonian occeig,pawmkrho,pawrhoij_alloc,pawrhoij_free,pawrhoij_init_unpacked pawrhoij_mpisum_unpacked,rf_transgrid_and_pack,sqnorm_v,symrhg,timab xmpi_sum
SOURCE
154 #if defined HAVE_CONFIG_H 155 #include "config.h" 156 #endif 157 158 #include "abi_common.h" 159 160 subroutine dfpt_vtorho(cg,cgq,cg1,cg1_active,cplex,cprj,cprjq,cprj1,dbl_nnsclo,& 161 & dim_eig2rf,doccde_rbz,docckqde,dtefield,dtfil,dtset,qphon,& 162 & edocc,eeig0,eigenq,eigen0,eigen1,ek0,ek1,eloc0,enl0,enl1,& 163 & fermie1,gh0c1_set,gh1c_set,gmet,gprimd,idir,indsy1,& 164 & ipert,irrzon1,istwfk_rbz,kg,kg1,kpt_rbz,mband,& 165 & mkmem,mkqmem,mk1mem,mpi_enreg,mpw,mpw1,my_natom,& 166 & natom,nband_rbz,ncpgr,nfftf,nhat1,nkpt_rbz,npwarr,npwar1,nres2,nspden,& 167 & nsppol,nsym1,ntypat,nvresid1,occkq,occ_rbz,optres,& 168 & paw_ij,paw_ij1,pawang,pawang1,pawfgr,pawfgrtab,pawrhoij,pawrhoij1,pawtab,& 169 & phnons1,ph1d,prtvol,psps,pwindall,qmat,resid,residm,rhog1,rhor1,rmet,rprimd,symaf1,symrc1,symrl1,ucvol,& 170 & usecprj,useylmgr1,ddk_f,vtrial,vtrial1,wtk_rbz,xred,ylm,ylm1,ylmgr1,cg1_out) 171 172 use defs_basis 173 use defs_datatypes 174 use defs_abitypes 175 use m_profiling_abi 176 use m_xmpi 177 use m_errors 178 use m_efield 179 use m_hamiltonian 180 use m_wfk 181 use m_cgtools 182 183 use m_occ, only : occeig 184 use m_hdr, only : hdr_skip, hdr_io 185 use m_pawang, only : pawang_type 186 use m_pawtab, only : pawtab_type 187 use m_paw_ij, only : paw_ij_type 188 use m_pawfgrtab,only : pawfgrtab_type 189 use m_pawrhoij, only : pawrhoij_type, pawrhoij_alloc, pawrhoij_free, & 190 & pawrhoij_init_unpacked, pawrhoij_free_unpacked, & 191 & pawrhoij_mpisum_unpacked 192 use m_pawcprj, only : pawcprj_type, pawcprj_alloc, pawcprj_free, pawcprj_get 193 use m_pawfgr, only : pawfgr_type 194 195 !This section has been created automatically by the script Abilint (TD). 196 !Do not modify the following lines by hand. 197 #undef ABI_FUNC 198 #define ABI_FUNC 'dfpt_vtorho' 199 use interfaces_18_timing 200 use interfaces_32_util 201 use interfaces_53_ffts 202 use interfaces_65_paw 203 use interfaces_66_wfs 204 use interfaces_67_common 205 use interfaces_72_response, except_this_one => dfpt_vtorho 206 !End of the abilint section 207 208 implicit none 209 210 !Arguments ------------------------------- 211 !scalars 212 integer,intent(in) :: cplex,dbl_nnsclo,dim_eig2rf,idir,ipert,mband,mk1mem,mkmem 213 integer,intent(in) :: mkqmem,mpw,mpw1,my_natom,natom,ncpgr,nfftf,nkpt_rbz,nspden 214 integer,intent(in) :: nsppol,nsym1,ntypat,optres,prtvol,usecprj,useylmgr1 215 integer,optional,intent(in) :: cg1_out 216 real(dp),intent(in) :: fermie1,ucvol 217 real(dp),intent(out) :: edocc,eeig0,ek0,ek1,eloc0,enl0,enl1,nres2,residm 218 type(MPI_type),intent(in) :: mpi_enreg 219 type(datafiles_type),intent(in) :: dtfil 220 type(dataset_type),intent(in) :: dtset 221 type(efield_type),intent(in) :: dtefield 222 type(pawang_type),intent(in) :: pawang,pawang1 223 type(pawfgr_type),intent(in) :: pawfgr 224 type(pseudopotential_type),intent(in) :: psps 225 226 !arrays 227 integer,intent(in) :: indsy1(4,nsym1,natom) 228 ! nfft**(1-1/nsym1) is 1 if nsym1==1, and nfft otherwise 229 integer,intent(in) :: irrzon1(dtset%nfft**(1-1/nsym1),2,(nspden/nsppol)-3*(nspden/4)) 230 integer,intent(in) :: istwfk_rbz(nkpt_rbz) 231 integer,intent(in) :: kg(3,mpw*mkmem),kg1(3,mpw1*mk1mem) 232 integer,intent(in) :: nband_rbz(nkpt_rbz*nsppol),npwar1(nkpt_rbz,2) 233 integer,intent(in) :: npwarr(nkpt_rbz,2),symaf1(nsym1),symrc1(3,3,nsym1),symrl1(3,3,nsym1) 234 real(dp),intent(in) :: qphon(3) 235 real(dp),intent(in) :: cg(2,mpw*dtset%nspinor*mband*mkmem*nsppol) 236 real(dp),intent(inout) :: cg1(2,mpw1*dtset%nspinor*mband*mk1mem*nsppol) 237 real(dp),intent(inout):: cg1_active(2,mpw1*dtset%nspinor*mband*mk1mem*nsppol*dim_eig2rf) 238 real(dp),intent(out) :: gh1c_set(2,mpw1*dtset%nspinor*mband*mk1mem*nsppol*dim_eig2rf) 239 real(dp),intent(out) :: gh0c1_set(2,mpw1*dtset%nspinor*mband*mk1mem*nsppol*dim_eig2rf) 240 real(dp),intent(in) :: cgq(2,mpw1*dtset%nspinor*mband*mkqmem*nsppol) 241 real(dp),intent(in) :: doccde_rbz(mband*nkpt_rbz*nsppol) 242 real(dp),intent(in) :: docckqde(mband*nkpt_rbz*nsppol) 243 real(dp),intent(in) :: eigen0(mband*nkpt_rbz*nsppol) 244 real(dp),intent(out) :: eigen1(2*mband*mband*nkpt_rbz*nsppol) 245 real(dp),intent(in) :: eigenq(mband*nkpt_rbz*nsppol),gmet(3,3),gprimd(3,3) 246 real(dp),intent(in) :: kpt_rbz(3,nkpt_rbz),occ_rbz(mband*nkpt_rbz*nsppol) 247 real(dp),intent(in) :: occkq(mband*nkpt_rbz*nsppol),ph1d(2,3*(2*dtset%mgfft+1)*natom) 248 real(dp),intent(in) :: phnons1(2,dtset%nfft**(1-1/nsym1),(nspden/nsppol)-3*(nspden/4)) 249 real(dp), intent(out) :: nhat1(cplex*nfftf,dtset%nspden*psps%usepaw) 250 real(dp),intent(out) :: resid(mband*nkpt_rbz*nsppol),rhog1(2,nfftf) 251 real(dp),intent(inout) :: nvresid1(cplex*nfftf,nspden),rhor1(cplex*nfftf,nspden) 252 real(dp),intent(in) :: rmet(3,3),rprimd(3,3) 253 real(dp),intent(in),target :: vtrial(nfftf,nspden) 254 real(dp),intent(inout),target :: vtrial1(cplex*nfftf,nspden) 255 real(dp),intent(in) :: wtk_rbz(nkpt_rbz),xred(3,natom) 256 real(dp),intent(in) :: ylm(mpw*mkmem,psps%mpsang*psps%mpsang*psps%useylm) 257 real(dp),intent(in) :: ylm1(mpw1*mk1mem,psps%mpsang*psps%mpsang*psps%useylm) 258 real(dp),intent(in) :: ylmgr1(mpw1*mk1mem,3+6*((ipert-natom)/10),psps%mpsang*psps%mpsang*psps%useylm*useylmgr1) 259 integer,intent(in) :: pwindall(max(mpw,mpw1)*mkmem,8,3) 260 real(dp),intent(in) :: qmat(2,dtefield%mband_occ,dtefield%mband_occ,nkpt_rbz,2,3) 261 type(pawcprj_type),intent(in) :: cprj (natom,dtset%nspinor*mband*mkmem *nsppol*usecprj) 262 type(pawcprj_type),intent(in) :: cprjq(natom,dtset%nspinor*mband*mkqmem*nsppol*usecprj) 263 type(pawcprj_type),intent(inout) :: cprj1(natom,dtset%nspinor*mband*mk1mem*nsppol*usecprj) 264 type(paw_ij_type),intent(in) :: paw_ij(my_natom*psps%usepaw),paw_ij1(my_natom*psps%usepaw) 265 type(pawfgrtab_type),intent(inout) :: pawfgrtab(my_natom*psps%usepaw) 266 type(pawrhoij_type),intent(in) :: pawrhoij(my_natom*psps%usepaw) 267 type(pawrhoij_type),target,intent(inout) :: pawrhoij1(my_natom*psps%usepaw) 268 type(pawtab_type), intent(in) :: pawtab(ntypat*psps%usepaw) 269 type(wfk_t),intent(inout) :: ddk_f(4) 270 271 !Local variables------------------------------- 272 !scalars 273 integer,parameter :: level=13 274 integer :: bd2tot_index,bdtot_index,buffer_size,counter,cplex_rhoij 275 integer :: iband,nlines_done,ibdkpt,ibg,ibg1,ibgq,icg,icg1,icgq,ierr 276 integer :: ii,ikg,ikg1,ikpt,ilm,index1,ispden,iscf_mod,isppol,istwf_k 277 integer :: mbd2kpsp,mbdkpsp,mcgq,mcgq_disk,mcprjq 278 integer :: mcprjq_disk,me,n1,n2,n3,n4,n5,n6,nband_k,nband_kq,nkpg,nkpg1 279 integer :: nnsclo_now,npw1_k,npw_k,nspden_rhoij,spaceworld,test_dot 280 logical :: paral_atom,qne0 281 real(dp) :: arg,wtk_k 282 type(gs_hamiltonian_type) :: gs_hamkq 283 type(rf_hamiltonian_type) :: rf_hamkq,rf_hamk_dir2 284 !arrays 285 integer,allocatable :: kg1_k(:,:),kg_k(:,:) 286 integer, pointer :: my_atmtab(:) 287 real(dp) :: kpoint(3),kpq(3) 288 real(dp) :: tsec(2) 289 real(dp),allocatable :: buffer1(:) 290 real(dp),allocatable :: ddkinpw(:),dkinpw(:),dkinpw2(:) 291 real(dp),allocatable :: doccde_k(:),doccde_kq(:) 292 real(dp),allocatable :: edocc_k(:),eeig0_k(:),eig0_k(:),eig0_kq(:),eig1_k(:) 293 real(dp),allocatable :: ek0_k(:),ek1_k(:),eloc0_k(:),enl0_k(:),enl1_k(:) 294 real(dp),allocatable :: ffnl1(:,:,:,:),ffnlk(:,:,:,:) 295 real(dp),allocatable :: grad_berry(:,:,:),kinpw1(:),kpg1_k(:,:) 296 real(dp),allocatable :: kpg_k(:,:),occ_k(:),occ_kq(:) 297 real(dp),allocatable :: ph3d(:,:,:),ph3d1(:,:,:),resid_k(:) 298 real(dp),allocatable :: rho1wfg(:,:),rho1wfr(:,:),rhoaug1(:,:,:,:),rocceig(:,:) 299 real(dp),allocatable :: vlocal(:,:,:,:),vlocal1(:,:,:,:) 300 real(dp),allocatable :: ylm1_k(:,:),ylm_k(:,:),ylmgr1_k(:,:,:) 301 type(pawrhoij_type),pointer :: pawrhoij1_unsym(:) 302 303 ! ********************************************************************* 304 305 DBG_ENTER('COLL') 306 307 !Keep track of total time spent in this routine 308 call timab(121,1,tsec) 309 call timab(124,1,tsec) 310 311 !Retrieve parallelism data 312 spaceworld=mpi_enreg%comm_cell 313 me=mpi_enreg%me_kpt 314 paral_atom=(my_natom/=natom) 315 my_atmtab=>mpi_enreg%my_atmtab 316 317 if (dtset%berryopt== 4.or.dtset%berryopt== 6.or.dtset%berryopt== 7.or.& 318 & dtset%berryopt==14.or.dtset%berryopt==16.or.dtset%berryopt==17) then 319 ABI_ALLOCATE(grad_berry,(2,mpw1,dtefield%mband_occ)) 320 else 321 ABI_ALLOCATE(grad_berry,(0,0,0)) 322 end if 323 324 !Test size of FFT grids (1 grid in norm-conserving, 2 grids in PAW) 325 if ((psps%usepaw==1.and.pawfgr%nfft/=nfftf).or.(psps%usepaw==0.and.dtset%nfft/=nfftf)) then 326 MSG_BUG('wrong values for nfft, nfftf!') 327 end if 328 329 !The value of iscf must be modified if ddk perturbation, see dfpt_looppert.f 330 iscf_mod=dtset%iscf;if(ipert==natom+1.or.ipert==natom+10.or.ipert==natom+11) iscf_mod=-3 331 332 edocc=zero ; eeig0=zero ; ek0=zero ; ek1=zero 333 eloc0=zero ; enl0=zero ; enl1=zero 334 bdtot_index=0 335 bd2tot_index=0 336 ibg=0;icg=0 337 ibgq=0;icgq=0 338 ibg1=0;icg1=0 339 mbdkpsp=mband*nkpt_rbz*nsppol 340 mbd2kpsp=2*mband**2*nkpt_rbz*nsppol 341 342 n1=dtset%ngfft(1); n2=dtset%ngfft(2); n3=dtset%ngfft(3) 343 n4=dtset%ngfft(4); n5=dtset%ngfft(5); n6=dtset%ngfft(6) 344 qne0=(qphon(1)**2+qphon(2)**2+qphon(3)**2>=tol14) 345 346 !Initialize PW 1st-order density if needed 347 !Also store old rho1 in case of density mixing 348 if (iscf_mod>0) then 349 if (optres==1) nvresid1=rhor1 350 if (psps%usepaw==0) then 351 rhor1(:,:)=zero 352 else 353 ABI_ALLOCATE(rho1wfr,(cplex*dtset%nfft,dtset%nspden)) 354 ABI_ALLOCATE(rho1wfg,(2,dtset%nfft)) 355 rho1wfr(:,:)=zero 356 end if 357 end if 358 359 !Set max number of non-self-consistent loops nnsclo_now for use in dfpt_vtowfk 360 if(iscf_mod<=0 .and. iscf_mod/=-3)then 361 nnsclo_now=dtset%nstep 362 else 363 if(dtset%nnsclo>0)then 364 nnsclo_now=dtset%nnsclo 365 else 366 nnsclo_now=1 367 end if 368 if(dbl_nnsclo==1) nnsclo_now=nnsclo_now*2 369 end if 370 371 !Prepare GS k+q wf 372 mcgq=mpw1*dtset%nspinor*mband*mkqmem*nsppol;mcgq_disk=0 373 374 !Prepare RF PAW files 375 if (psps%usepaw==1) then 376 mcprjq=dtset%nspinor*mband*mkqmem*nsppol*usecprj;mcprjq_disk=0 377 else 378 mcprjq=0;mcprjq_disk=0 379 end if 380 381 !Initialisation of the wfdot file in case of electric field (or 2nd order Sternheimer equation) 382 test_dot=0 383 if (ipert==natom+2.and.sum((dtset%qptn(1:3))**2 )<=tol7.and.& 384 & (dtset%berryopt/= 4.and.dtset%berryopt/= 6.and.dtset%berryopt/= 7.and.& 385 & dtset%berryopt/=14.and.dtset%berryopt/=16.and.dtset%berryopt/=17).or.& 386 & (ipert==natom+10.or.ipert==natom+11)) then 387 test_dot=1 388 end if 389 390 !==== Initialize most of the Hamiltonian (and derivative) ==== 391 !1) Allocate all arrays and initialize quantities that do not depend on k and spin. 392 !2) Perform the setup needed for the non-local factors: 393 !* Norm-conserving: Constant kleimann-Bylander energies are copied from psps to gs_hamk. 394 !* PAW: Initialize the overlap coefficients and allocate the Dij coefficients. 395 396 call init_hamiltonian(gs_hamkq,psps,pawtab,dtset%nspinor,nsppol,nspden,natom,& 397 & dtset%typat,xred,dtset%nfft,dtset%mgfft,dtset%ngfft,rprimd,dtset%nloalg,& 398 & paw_ij=paw_ij,comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab,mpi_spintab=mpi_enreg%my_isppoltab,& 399 & usecprj=usecprj,ph1d=ph1d,nucdipmom=dtset%nucdipmom,use_gpu_cuda=dtset%use_gpu_cuda) 400 401 call init_rf_hamiltonian(cplex,gs_hamkq,ipert,rf_hamkq,has_e1kbsc=.true.,paw_ij1=paw_ij1,& 402 & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab,mpi_spintab=mpi_enreg%my_isppoltab) 403 if ((ipert==natom+10.and.idir>3).or.ipert==natom+11) then 404 call init_rf_hamiltonian(cplex,gs_hamkq,ipert,rf_hamk_dir2,has_e1kbsc=.true.,paw_ij1=paw_ij1,& 405 & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab,mpi_spintab=mpi_enreg%my_isppoltab) 406 end if 407 408 !PAW:allocate memory for non-symetrized 1st-order occupancies matrix (pawrhoij1) 409 pawrhoij1_unsym => pawrhoij1 410 if (psps%usepaw==1.and.iscf_mod>0) then 411 if (paral_atom) then 412 ABI_DATATYPE_ALLOCATE(pawrhoij1_unsym,(natom)) 413 cplex_rhoij=max(cplex,dtset%pawcpxocc);nspden_rhoij=dtset%nspden 414 call pawrhoij_alloc(pawrhoij1_unsym,cplex_rhoij,nspden_rhoij,dtset%nspinor,& 415 & dtset%nsppol,dtset%typat,pawtab=pawtab,use_rhoijp=0,use_rhoij_=1) 416 else 417 pawrhoij1_unsym => pawrhoij1 418 call pawrhoij_init_unpacked(pawrhoij1_unsym) 419 end if 420 end if 421 422 ABI_ALLOCATE(rhoaug1,(cplex*n4,n5,n6,gs_hamkq%nvloc)) 423 ABI_ALLOCATE(vlocal,(n4,n5,n6,gs_hamkq%nvloc)) 424 ABI_ALLOCATE(vlocal1,(cplex*n4,n5,n6,gs_hamkq%nvloc)) 425 426 nlines_done = 0 427 428 !LOOP OVER SPINS 429 do isppol=1,nsppol 430 431 ! Rewind kpgsph data file if needed: 432 ikg=0;ikg1=0 433 434 ! Set up local potential vlocal1 with proper dimensioning, from vtrial1 435 ! Same thing for vlocal from vtrial Also take into account the spin. 436 437 call rf_transgrid_and_pack(isppol,nspden,psps%usepaw,cplex,nfftf,dtset%nfft,dtset%ngfft,& 438 & gs_hamkq%nvloc,pawfgr,mpi_enreg,vtrial,vtrial1,vlocal,vlocal1) 439 440 ! Continue to initialize the Hamiltonian 441 call load_spin_hamiltonian(gs_hamkq,isppol,vlocal=vlocal,with_nonlocal=.true.) 442 call load_spin_rf_hamiltonian(rf_hamkq,isppol,vlocal1=vlocal1,with_nonlocal=.true.) 443 if ((ipert==natom+10.and.idir>3).or.ipert==natom+11) then 444 call load_spin_rf_hamiltonian(rf_hamk_dir2,isppol,with_nonlocal=.true.) 445 if (ipert==natom+11) then ! load vlocal1 446 call load_spin_rf_hamiltonian(rf_hamk_dir2,isppol,vlocal1=vlocal1) 447 end if 448 end if 449 450 if (ipert==natom+5) then !SPr deb, in case of magnetic field perturbation, no non-local 451 call load_spin_rf_hamiltonian(rf_hamkq,isppol,vlocal1=vlocal1) 452 end if 453 454 ! Nullify contribution to 1st-order density from this k-point 455 rhoaug1(:,:,:,:)=zero 456 457 call timab(125,1,tsec) 458 459 !====================================================================== 460 !============== BIG FAT K POINT LOOP ================================ 461 !====================================================================== 462 463 do ikpt=1,nkpt_rbz 464 counter=100*ikpt+isppol 465 466 nband_k = nband_rbz(ikpt+(isppol-1)*nkpt_rbz) 467 istwf_k = istwfk_rbz(ikpt) 468 npw_k = npwarr(ikpt,1) 469 npw1_k = npwar1(ikpt,1) 470 wtk_k = wtk_rbz(ikpt) 471 472 if(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,isppol,me)) then 473 eigen1(1+bd2tot_index : 2*nband_k**2+bd2tot_index) = zero 474 resid(1+bdtot_index : nband_k+bdtot_index) = zero 475 bdtot_index=bdtot_index+nband_k 476 bd2tot_index=bd2tot_index+2*nband_k**2 477 478 cycle ! Skip the rest of the k-point loop 479 end if 480 481 kpoint(:)=kpt_rbz(:,ikpt) 482 kpq(:)=kpoint(:);if (ipert<natom+3.or.ipert==natom+5) kpq(:)=kpq(:)+qphon(1:3) 483 ABI_ALLOCATE(kg_k,(3,npw_k)) 484 ABI_ALLOCATE(kg1_k,(3,npw1_k)) 485 ABI_ALLOCATE(ylm_k,(npw_k,psps%mpsang*psps%mpsang*psps%useylm)) 486 ABI_ALLOCATE(ylm1_k,(npw1_k,psps%mpsang*psps%mpsang*psps%useylm)) 487 ABI_ALLOCATE(ylmgr1_k,(npw1_k,3+6*((ipert-natom)/10),psps%mpsang*psps%mpsang*psps%useylm*useylmgr1)) 488 ABI_ALLOCATE(doccde_k,(nband_k)) 489 ABI_ALLOCATE(doccde_kq,(nband_k)) 490 ABI_ALLOCATE(eig0_k,(nband_k)) 491 ABI_ALLOCATE(eig0_kq,(nband_k)) 492 ABI_ALLOCATE(eig1_k,(2*nband_k**2)) 493 ABI_ALLOCATE(edocc_k,(nband_k)) 494 ABI_ALLOCATE(eeig0_k,(nband_k)) 495 ABI_ALLOCATE(ek0_k,(nband_k)) 496 ABI_ALLOCATE(ek1_k,(nband_k)) 497 ABI_ALLOCATE(eloc0_k,(nband_k)) 498 ABI_ALLOCATE(occ_k,(nband_k)) 499 ABI_ALLOCATE(occ_kq,(nband_k)) 500 ABI_ALLOCATE(resid_k,(nband_k)) 501 ABI_ALLOCATE(rocceig,(nband_k,nband_k)) 502 ABI_ALLOCATE(enl0_k,(nband_k)) 503 ABI_ALLOCATE(enl1_k,(nband_k)) 504 505 eig1_k(:)=zero 506 eig0_k(:)=eigen0(1+bdtot_index:nband_k+bdtot_index) 507 eig0_kq(:)=eigenq(1+bdtot_index:nband_k+bdtot_index) 508 edocc_k(:)=zero 509 eeig0_k(:)=zero ; ek0_k(:)=zero ; ek1_k(:)=zero 510 eloc0_k(:)=zero ; enl0_k(:)=zero ; enl1_k(:)=zero 511 occ_k(:)=occ_rbz(1+bdtot_index:nband_k+bdtot_index) 512 occ_kq(:)=occkq(1+bdtot_index:nband_k+bdtot_index) 513 doccde_k(:)=doccde_rbz(1+bdtot_index:nband_k+bdtot_index) 514 doccde_kq(:)=docckqde(1+bdtot_index:nband_k+bdtot_index) 515 resid_k(:)=zero 516 517 ! For each pair of active bands (m,n), generates the ratios 518 ! rocceig(m,n)=(occ_kq(m)-occ_k(n))/(eig0_kq(m)-eig0_k(n)) 519 ! and decide to which band to attribute it. 520 call occeig(doccde_k,doccde_kq,eig0_k,eig0_kq,nband_k,dtset%occopt,occ_k,occ_kq,rocceig) 521 522 ! These arrays are not needed anymore. 523 ABI_DEALLOCATE(doccde_k) 524 ABI_DEALLOCATE(doccde_kq) 525 ABI_DEALLOCATE(occ_kq) 526 527 kg_k(:,1:npw_k)=kg(:,1+ikg:npw_k+ikg) 528 if (psps%useylm==1) then 529 do ilm=1,psps%mpsang*psps%mpsang 530 ylm_k(1:npw_k,ilm)=ylm(1+ikg:npw_k+ikg,ilm) 531 end do 532 end if 533 534 ! Get (k+q+G) wave vectors and associated spherical harmonics 535 kg1_k(:,1:npw1_k)=kg1(:,1+ikg1:npw1_k+ikg1) 536 if (psps%useylm==1) then 537 do ilm=1,psps%mpsang*psps%mpsang 538 ylm1_k(1:npw1_k,ilm)=ylm1(1+ikg1:npw1_k+ikg1,ilm) 539 end do 540 if (useylmgr1==1) then 541 do ilm=1,psps%mpsang*psps%mpsang 542 do ii=1,3+6*((ipert-natom)/10) 543 ylmgr1_k(1:npw1_k,ii,ilm)=ylmgr1(1+ikg1:npw1_k+ikg1,ii,ilm) 544 end do 545 end do 546 end if 547 end if 548 549 ! Set up the ground-state Hamiltonian, and some parts of the 1st-order Hamiltonian 550 call getgh1c_setup(gs_hamkq,rf_hamkq,dtset,psps,& ! In 551 kpoint,kpq,idir,ipert,natom,rmet,gprimd,gmet,istwf_k,& ! In 552 npw_k,npw1_k,useylmgr1,kg_k,ylm_k,kg1_k,ylm1_k,ylmgr1_k,& ! In 553 dkinpw,nkpg,nkpg1,kpg_k,kpg1_k,kinpw1,ffnlk,ffnl1,ph3d,ph3d1,& ! Out 554 ddkinpw=ddkinpw,dkinpw2=dkinpw2,rf_hamk_dir2=rf_hamk_dir2) ! Out 555 556 ! Compute the gradient of the Berry-phase term 557 if (dtset%berryopt== 4.or.dtset%berryopt== 6.or.dtset%berryopt== 7.or.& 558 & dtset%berryopt==14.or.dtset%berryopt==16.or.dtset%berryopt==17) then 559 if (ipert<=natom) then 560 ! phonon perturbation 561 call dfptff_gradberry(cg,cg1,dtefield,grad_berry,ikpt,isppol,mband,mpw,mpw1,mkmem,mk1mem,nkpt_rbz,& 562 & npwarr,npwar1,dtset%nspinor,nsppol,qmat,pwindall) 563 else 564 ! electric field perturbation 565 call dfptff_gbefd(cg,cg1,dtefield,grad_berry,idir,ikpt,isppol,mband,mpw,mpw1,mkmem,mk1mem,nkpt_rbz,& 566 & npwarr,npwar1,dtset%nspinor,& 567 & nsppol,qmat,pwindall,rprimd) 568 end if 569 end if 570 571 ! Free some memory before calling dfpt_vtowfk 572 ABI_DEALLOCATE(ylm_k) 573 ABI_DEALLOCATE(ylm1_k) 574 ABI_DEALLOCATE(ylmgr1_k) 575 576 ! Compute the eigenvalues, wavefunction, residuals, 577 ! contributions to kinetic energy, nonlocal energy, forces, 578 ! and update of 1st-order density to this k-point and this spin polarization. 579 nband_kq = nband_k !Note that the calculation only works for same number of bandes on all K points. 580 ! Note that dfpt_vtowfk is called with kpoint, while kpt is used inside vtowfk3 581 call dfpt_vtowfk(cg,cgq,cg1,cg1_active,cplex,cprj,cprjq,cprj1,dim_eig2rf,dtfil,& 582 & dtset,edocc_k,eeig0_k,eig0_k,eig0_kq,eig1_k,ek0_k,ek1_k,eloc0_k,enl0_k,enl1_k,fermie1,& 583 & gh0c1_set,gh1c_set,grad_berry,gs_hamkq,ibg,ibgq,ibg1,icg,icgq,icg1,idir,ikpt,ipert,isppol,& 584 & mband,mcgq,mcprjq,mkmem,mk1mem,mpi_enreg,mpw,mpw1,natom,nband_k,ncpgr,nnsclo_now,& 585 & npw_k,npw1_k,dtset%nspinor,nsppol,n4,n5,n6,occ_k,pawrhoij1_unsym,prtvol,psps,resid_k,& 586 & rf_hamkq,rf_hamk_dir2,rhoaug1,rocceig,ddk_f,wtk_k,nlines_done,cg1_out) 587 588 ! Free temporary storage 589 ABI_DEALLOCATE(kinpw1) 590 ABI_DEALLOCATE(kg_k) 591 ABI_DEALLOCATE(kg1_k) 592 ABI_DEALLOCATE(kpg_k) 593 ABI_DEALLOCATE(kpg1_k) 594 ABI_DEALLOCATE(dkinpw) 595 if (ipert==natom+10) then 596 ABI_DEALLOCATE(ddkinpw) 597 if (idir>3) then 598 ABI_DEALLOCATE(dkinpw2) 599 end if 600 end if 601 ABI_DEALLOCATE(ffnlk) 602 ABI_DEALLOCATE(ffnl1) 603 ABI_DEALLOCATE(eig0_k) 604 ABI_DEALLOCATE(eig0_kq) 605 ABI_DEALLOCATE(rocceig) 606 ABI_DEALLOCATE(ph3d) 607 if (allocated(ph3d1)) then 608 ABI_DEALLOCATE(ph3d1) 609 end if 610 611 ! Save eigenvalues (hartree), residuals (hartree**2) 612 eigen1 (1+bd2tot_index : 2*nband_k**2+bd2tot_index) = eig1_k(:) 613 resid (1+bdtot_index : nband_k+bdtot_index) = resid_k(:) 614 615 ! Accumulate sum over k points for nonlocal and kinetic energies, 616 ! also accumulate gradients of Enonlocal: 617 if (iscf_mod>0 .or. iscf_mod==-3 .or. iscf_mod==-2)then 618 do iband=1,nband_k 619 edocc=edocc+wtk_k*occ_k(iband)*edocc_k(iband) 620 eeig0=eeig0+wtk_k*occ_k(iband)*eeig0_k(iband) 621 ek0=ek0+wtk_k*occ_k(iband)*ek0_k(iband) 622 ek1=ek1+wtk_k*occ_k(iband)*ek1_k(iband) 623 eloc0=eloc0+wtk_k*occ_k(iband)*eloc0_k(iband) 624 enl0=enl0+wtk_k*occ_k(iband)*enl0_k(iband) 625 enl1=enl1+wtk_k*occ_k(iband)*enl1_k(iband) 626 end do 627 end if 628 629 ABI_DEALLOCATE(eig1_k) 630 ABI_DEALLOCATE(occ_k) 631 ABI_DEALLOCATE(resid_k) 632 ABI_DEALLOCATE(edocc_k) 633 ABI_DEALLOCATE(eeig0_k) 634 ABI_DEALLOCATE(ek0_k) 635 ABI_DEALLOCATE(ek1_k) 636 ABI_DEALLOCATE(eloc0_k) 637 ABI_DEALLOCATE(enl0_k) 638 ABI_DEALLOCATE(enl1_k) 639 640 ! Keep track of total number of bands (all k points so far, even for k points not treated by me) 641 bdtot_index=bdtot_index+nband_k 642 bd2tot_index=bd2tot_index+2*nband_k**2 643 644 ! Shift array memory 645 if (mkmem/=0) then 646 ibg=ibg+nband_k 647 icg=icg+npw_k*dtset%nspinor*nband_k 648 ikg=ikg+npw_k 649 end if 650 if (mkqmem/=0) then 651 ibgq=ibgq+dtset%nspinor*nband_k 652 icgq=icgq+npw1_k*dtset%nspinor*nband_k 653 end if 654 if (mk1mem/=0) then 655 ibg1=ibg1+nband_k 656 icg1=icg1+npw1_k*dtset%nspinor*nband_k 657 ikg1=ikg1+npw1_k 658 end if 659 660 end do 661 662 !====================================================================== 663 !================== END BIG K POINT LOOP ============================ 664 !====================================================================== 665 666 call timab(125,2,tsec) 667 668 ! Transfer density on augmented fft grid to normal fft grid in real space. Also take into account the spin. 669 ! FR EB for the non-collinear part see vtorho.F90 670 if(iscf_mod>0) then 671 if (psps%usepaw==0) then 672 call fftpac(isppol,mpi_enreg,nspden,cplex*n1,n2,n3,cplex*n4,n5,n6,dtset%ngfft,rhor1,rhoaug1(:,:,:,1),1) 673 if(nspden==4)then 674 do ispden=2,4 675 call fftpac(ispden,mpi_enreg,nspden,cplex*n1,n2,n3,cplex*n4,n5,n6,dtset%ngfft,rhor1,rhoaug1(:,:,:,ispden),1) 676 end do 677 end if 678 else 679 call fftpac(isppol,mpi_enreg,nspden,cplex*n1,n2,n3,cplex*n4,n5,n6,dtset%ngfft,rho1wfr,rhoaug1(:,:,:,1),1) 680 if(nspden==4)then 681 do ispden=2,4 682 call fftpac(ispden,mpi_enreg,nspden,cplex*n1,n2,n3,cplex*n4,n5,n6,dtset%ngfft,rho1wfr,rhoaug1(:,:,:,ispden),1) 683 end do 684 end if 685 end if 686 end if 687 688 end do ! End loop over spins 689 690 !More memory cleaning 691 call destroy_hamiltonian(gs_hamkq) 692 call destroy_rf_hamiltonian(rf_hamkq) 693 if ((ipert==natom+10.and.idir>3).or.ipert==natom+11) then 694 call destroy_rf_hamiltonian(rf_hamk_dir2) 695 end if 696 ABI_DEALLOCATE(rhoaug1) 697 ABI_DEALLOCATE(vlocal) 698 ABI_DEALLOCATE(vlocal1) 699 700 call timab(124,2,tsec) 701 702 !=== MPI communications ================== 703 if(xmpi_paral==1)then 704 call timab(129,1,tsec) 705 706 ! Compute buffer size 707 buffer_size=7+mbd2kpsp+mbdkpsp 708 if (iscf_mod>0) then 709 buffer_size=buffer_size+cplex*dtset%nfft*nspden 710 end if 711 ABI_ALLOCATE(buffer1,(buffer_size)) 712 713 ! Pack rhor1,edocc,eeig0,ek0,ek1,eloc0,enl0,enl1,eigen1,resid 714 if (iscf_mod>0) then 715 index1=cplex*dtset%nfft*nspden 716 if (psps%usepaw==0) then 717 buffer1(1:index1)=reshape(rhor1 ,(/index1/)) 718 else 719 buffer1(1:index1)=reshape(rho1wfr,(/index1/)) 720 end if 721 else 722 index1=0 723 end if 724 buffer1(index1+1)=edocc;buffer1(index1+2)=eeig0 725 buffer1(index1+3)=ek0 ;buffer1(index1+4)=ek1 726 buffer1(index1+5)=eloc0;buffer1(index1+6)=enl0 727 buffer1(index1+7)=enl1 728 index1=index1+7 729 bdtot_index=0;bd2tot_index=0 730 do isppol=1,nsppol 731 do ikpt=1,nkpt_rbz 732 nband_k=nband_rbz(ikpt+(isppol-1)*nkpt_rbz) 733 buffer1(index1+1:index1+2*nband_k**2) = eigen1(bd2tot_index+1:bd2tot_index+2*nband_k**2) 734 buffer1(index1+2*nband_k**2+1:index1+2*nband_k**2+nband_k)= resid(bdtot_index+1:bdtot_index+nband_k) 735 bdtot_index=bdtot_index+nband_k 736 bd2tot_index=bd2tot_index+2*nband_k**2 737 index1=index1+2*nband_k**2+nband_k 738 end do 739 end do 740 if(index1<buffer_size)buffer1(index1+1:buffer_size)=zero 741 742 ! Build sum of everything 743 call timab(48,1,tsec) 744 call xmpi_sum(buffer1,buffer_size,spaceworld,ierr) 745 call timab(48,2,tsec) 746 747 ! Unpack the final result 748 if(iscf_mod>0) then 749 index1=cplex*dtset%nfft*nspden 750 if (psps%usepaw==0) then 751 rhor1(:,:) =reshape(buffer1(1:index1),(/cplex*dtset%nfft,nspden/)) 752 else 753 rho1wfr(:,:)=reshape(buffer1(1:index1),(/cplex*dtset%nfft,nspden/)) 754 end if 755 else 756 index1=0 757 end if 758 759 edocc=buffer1(index1+1);eeig0=buffer1(index1+2) 760 ek0=buffer1(index1+3) ;ek1=buffer1(index1+4) 761 eloc0=buffer1(index1+5);enl0=buffer1(index1+6) 762 enl1=buffer1(index1+7) 763 index1=index1+7 764 bdtot_index=0;bd2tot_index=0 765 do isppol=1,nsppol 766 do ikpt=1,nkpt_rbz 767 nband_k=nband_rbz(ikpt+(isppol-1)*nkpt_rbz) 768 eigen1(bd2tot_index+1:bd2tot_index+2*nband_k**2) = buffer1(index1+1:index1+2*nband_k**2) 769 resid(bdtot_index+1:bdtot_index+nband_k)= buffer1(index1+2*nband_k**2+1:index1+2*nband_k**2+nband_k) 770 bdtot_index=bdtot_index+nband_k 771 bd2tot_index=bd2tot_index+2*nband_k**2 772 index1=index1+2*nband_k**2+nband_k 773 end do 774 end do 775 ABI_DEALLOCATE(buffer1) 776 777 ! Accumulate PAW occupancies 778 if (psps%usepaw==1.and.iscf_mod>0) then 779 call pawrhoij_mpisum_unpacked(pawrhoij1_unsym,spaceworld) 780 end if 781 782 call timab(129,2,tsec) 783 end if ! if kpt parallel 784 785 call timab(127,1,tsec) 786 787 !If needed, compute rhog1, and symmetrizes the density 788 if (iscf_mod > 0) then 789 790 ! In order to have the symrhg working in parallel on FFT coefficients, the size 791 ! of irzzon1 and phnons1 should be set to nfftot. Therefore, nsym\=1 does not work. 792 793 if(nspden==4) then 794 ! FR symrhg will manage correctly this rearrangement 795 rhor1(:,2)=rhor1(:,2)+(rhor1(:,1)+rhor1(:,4)) !(n+mx) 796 rhor1(:,3)=rhor1(:,3)+(rhor1(:,1)+rhor1(:,4)) !(n+my) 797 call timab(17,2,tsec) 798 end if 799 ! 800 if (psps%usepaw==0) then 801 call symrhg(cplex,gprimd,irrzon1,mpi_enreg,dtset%nfft,dtset%nfft,dtset%ngfft,& 802 & nspden,nsppol,nsym1,dtset%paral_kgb,phnons1,rhog1 ,rhor1 ,rprimd,symaf1,symrl1) 803 else 804 call symrhg(cplex,gprimd,irrzon1,mpi_enreg,dtset%nfft,dtset%nfft,dtset%ngfft,& 805 & nspden,nsppol,nsym1,dtset%paral_kgb,phnons1,rho1wfg,rho1wfr,rprimd,symaf1,symrl1) 806 end if 807 ! We now have both rho(r) and rho(G), symmetrized, and if nsppol=2 808 ! we also have the spin-up density, symmetrized, in rhor1(:,2). 809 end if 810 811 ABI_DEALLOCATE(grad_berry) 812 813 !Find largest residual over bands, k points, and spins except for nbdbuf highest bands 814 ibdkpt=1 815 residm=zero 816 do isppol=1,nsppol 817 do ikpt=1,nkpt_rbz 818 nband_k=nband_rbz(ikpt+(isppol-1)*nkpt_rbz) 819 nband_k=max(1,nband_k-dtset%nbdbuf) 820 residm=max(residm,maxval(resid(ibdkpt:ibdkpt+nband_k-1))) 821 ibdkpt=ibdkpt+nband_k 822 end do 823 end do 824 825 call timab(127,2,tsec) 826 827 if (iscf_mod>0) then 828 829 ! PAW: Build new 1st-order rhoij quantities then symetrize them 830 ! Compute and add the 1st-order compensation density to rho1wfr 831 ! to get the total 1st-order density 832 if (psps%usepaw==1) then 833 call pawmkrho(arg,cplex,gprimd,idir,indsy1,ipert,mpi_enreg,& 834 & my_natom,natom,nspden,nsym1,ntypat,dtset%paral_kgb,pawang,pawfgr,pawfgrtab,& 835 & dtset%pawprtvol,pawrhoij1,pawrhoij1_unsym,pawtab,dtset%qptn,rho1wfg,rho1wfr,& 836 & rhor1,rprimd,symaf1,symrc1,dtset%typat,ucvol,dtset%usewvl,xred,& 837 & pawang_sym=pawang1,pawnhat=nhat1,pawrhoij0=pawrhoij,rhog=rhog1) 838 ABI_DEALLOCATE(rho1wfr) 839 ABI_DEALLOCATE(rho1wfg) 840 if (paral_atom) then 841 call pawrhoij_free(pawrhoij1_unsym) 842 ABI_DATATYPE_DEALLOCATE(pawrhoij1_unsym) 843 end if 844 end if 845 846 ! Compute density residual (if required) and its squared norm 847 if (optres==1) then 848 nvresid1=rhor1-nvresid1 849 call sqnorm_v(1,nfftf,nres2,dtset%nspden,optres,nvresid1) 850 end if 851 end if ! iscf>0 852 853 call timab(121,2,tsec) 854 855 DBG_EXIT('COLL') 856 857 end subroutine dfpt_vtorho