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
INPUTS
cg(2,mpw*nspinor*mband_mem*mkmem*nsppol)=planewave coefficients of wavefunctions cgq(2,mpw1*nspinor*mband_mem*mkqmem*nsppol)=pw coefficients of GS wavefunctions at k+q. cg1(2,mpw1*nspinor*mband_mem*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 mband_mem=maximum number of bands on this cpu 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 tnons1(3,nsym1)=non-symmorphic translations 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 vectornd(with_vectornd*nfftf,nspden,3)=nuclear dipole moment vector potential vtrial(nfftf,nspden)=GS Vtrial(r). vtrial1(cplex*nfftf,nspden)=INPUT RF Vtrial(r). vxctau(nfftf,nspden,4*usekden)=derivative of e_xc wrt kin energy density, for mGGA with_vectornd = 1 if vectornd allocated 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_mem*mk1mem*nsppol)=updated wavefunctions, orthogonalized to the occupied states
cg1_active(2,mpw1*nspinor*mband_mem*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
end0=0th-order nuclear dipole energy part of 2nd-order total energy.
end1=1st-order nuclear dipole energy 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.
evxctau0=0th-order energy from vxctau
evxctau1=1st-order energy from vxctau
gh1c_set(2,mpw1*nspinor*mband_mem*mk1mem*nsppol*dim_eig2rf)= set of <G|H^{(1)}|nK>
gh0c1_set(2,mpw1*nspinor*mband_mem*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_mem*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.
SOURCE
213 subroutine dfpt_vtorho(cg,cgq,cg1,cg1_active,cplex,cprj,cprjq,cprj1,dbl_nnsclo,& 214 & dim_eig2rf,doccde_rbz,docckqde,dtefield,dtfil,dtset,qphon,& 215 & edocc,eeig0,eigenq,eigen0,eigen1,ek0,ek1,eloc0,end0,end1,enl0,enl1,& 216 & evxctau0,evxctau1,fermie1,gh0c1_set,gh1c_set,gmet,gprimd,idir,indsy1,& 217 & ipert,irrzon1,istwfk_rbz,kg,kg1,kpt_rbz,mband,mband_mem,& 218 & mkmem,mkqmem,mk1mem,mpi_enreg,mpw,mpw1,my_natom,& 219 & natom,nband_rbz,ncpgr,nfftf,nhat1,nkpt_rbz,npwarr,npwar1,nres2,nspden,& 220 & nsppol,nsym1,ntypat,nvresid1,occkq,occ_rbz,optres,& 221 & paw_ij,paw_ij1,pawang,pawang1,pawfgr,pawfgrtab,pawrhoij,pawrhoij1,pawtab,& 222 & phnons1,ph1d,prtvol,psps,pwindall,qmat,resid,residm,rhog1,rhor1,rmet,rprimd,symaf1,symrc1,symrl1,tnons1,ucvol,& 223 & usecprj,useylmgr1,ddk_f,vectornd,vtrial,vtrial1,vxctau,with_vectornd,wtk_rbz,xred,ylm,ylm1,ylmgr1,cg1_out) 224 225 !Arguments ------------------------------- 226 !scalars 227 integer,intent(in) :: cplex,dbl_nnsclo,dim_eig2rf,idir,ipert,mband,mk1mem,mkmem 228 integer,intent(in) :: mband_mem 229 integer,intent(in) :: mkqmem,mpw,mpw1,my_natom,natom,ncpgr,nfftf,nkpt_rbz,nspden 230 integer,intent(in) :: nsppol,nsym1,ntypat,optres,prtvol,usecprj,useylmgr1,with_vectornd 231 integer,optional,intent(in) :: cg1_out 232 real(dp),intent(in) :: fermie1,ucvol 233 real(dp),intent(out) :: edocc,eeig0,ek0,ek1,eloc0,end0,end1,enl0,enl1,evxctau0,evxctau1,nres2,residm 234 type(MPI_type),intent(in) :: mpi_enreg 235 type(datafiles_type),intent(in) :: dtfil 236 type(dataset_type),intent(in) :: dtset 237 type(efield_type),intent(in) :: dtefield 238 type(pawang_type),intent(in) :: pawang,pawang1 239 type(pawfgr_type),intent(in) :: pawfgr 240 type(pseudopotential_type),intent(in) :: psps 241 242 !arrays 243 integer,intent(in) :: indsy1(4,nsym1,natom) 244 ! nfft**(1-1/nsym1) is 1 if nsym1==1, and nfft otherwise 245 integer,intent(in) :: irrzon1(dtset%nfft**(1-1/nsym1),2,(nspden/nsppol)-3*(nspden/4)) 246 integer,intent(in) :: istwfk_rbz(nkpt_rbz) 247 integer,intent(in) :: kg(3,mpw*mkmem),kg1(3,mpw1*mk1mem) 248 integer,intent(in) :: nband_rbz(nkpt_rbz*nsppol),npwar1(nkpt_rbz,2) 249 integer,intent(in) :: npwarr(nkpt_rbz,2),symaf1(nsym1),symrc1(3,3,nsym1),symrl1(3,3,nsym1) 250 real(dp),intent(in) :: qphon(3) 251 real(dp),intent(in) :: cg(2,mpw*dtset%nspinor*mband_mem*mkmem*nsppol) 252 real(dp),intent(inout) :: cg1(2,mpw1*dtset%nspinor*mband_mem*mk1mem*nsppol) 253 real(dp),intent(inout):: cg1_active(2,mpw1*dtset%nspinor*mband_mem*mk1mem*nsppol*dim_eig2rf) 254 real(dp),intent(out) :: gh1c_set(2,mpw1*dtset%nspinor*mband_mem*mk1mem*nsppol*dim_eig2rf) 255 real(dp),intent(out) :: gh0c1_set(2,mpw1*dtset%nspinor*mband_mem*mk1mem*nsppol*dim_eig2rf) 256 real(dp),intent(in) :: cgq(2,mpw1*dtset%nspinor*mband_mem*mkqmem*nsppol) 257 real(dp),intent(in) :: doccde_rbz(mband*nkpt_rbz*nsppol) 258 real(dp),intent(in) :: docckqde(mband*nkpt_rbz*nsppol) 259 real(dp),intent(in) :: eigen0(mband*nkpt_rbz*nsppol) 260 real(dp),intent(out) :: eigen1(2*mband*mband*nkpt_rbz*nsppol) 261 real(dp),intent(in) :: eigenq(mband*nkpt_rbz*nsppol),gmet(3,3),gprimd(3,3) 262 real(dp),intent(in) :: kpt_rbz(3,nkpt_rbz),occ_rbz(mband*nkpt_rbz*nsppol) 263 real(dp),intent(in) :: occkq(mband*nkpt_rbz*nsppol),ph1d(2,3*(2*dtset%mgfft+1)*natom) 264 real(dp),intent(in) :: phnons1(2,dtset%nfft**(1-1/nsym1),(nspden/nsppol)-3*(nspden/4)) 265 real(dp), intent(out) :: nhat1(cplex*nfftf,dtset%nspden*psps%usepaw) 266 real(dp),intent(out) :: resid(mband*nkpt_rbz*nsppol),rhog1(2,nfftf) 267 real(dp),intent(inout) :: nvresid1(cplex*nfftf,nspden),rhor1(cplex*nfftf,nspden) 268 real(dp),intent(inout) :: vectornd(with_vectornd*nfftf,dtset%nspden,3) 269 real(dp),intent(in) :: rmet(3,3),rprimd(3,3) 270 real(dp),intent(in) :: tnons1(3,nsym1) 271 real(dp),intent(in),target :: vtrial(nfftf,nspden) 272 real(dp),intent(inout),target :: vtrial1(cplex*nfftf,nspden) 273 real(dp),intent(inout) :: vxctau(nfftf,dtset%nspden,4*dtset%usekden) 274 real(dp),intent(in) :: wtk_rbz(nkpt_rbz),xred(3,natom) 275 real(dp),intent(in) :: ylm(mpw*mkmem,psps%mpsang*psps%mpsang*psps%useylm) 276 real(dp),intent(in) :: ylm1(mpw1*mk1mem,psps%mpsang*psps%mpsang*psps%useylm) 277 real(dp),intent(in) :: ylmgr1(mpw1*mk1mem,3+6*((ipert-natom)/10),psps%mpsang*psps%mpsang*psps%useylm*useylmgr1) 278 integer,intent(in) :: pwindall(max(mpw,mpw1)*mkmem,8,3) 279 real(dp),intent(in) :: qmat(2,dtefield%mband_occ,dtefield%mband_occ,nkpt_rbz,2,3) 280 type(pawcprj_type),intent(in) :: cprj (natom,dtset%nspinor*mband_mem*mkmem *nsppol*usecprj) 281 type(pawcprj_type),intent(in) :: cprjq(natom,dtset%nspinor*mband_mem*mkqmem*nsppol*usecprj) 282 type(pawcprj_type),intent(inout) :: cprj1(natom,dtset%nspinor*mband_mem*mk1mem*nsppol*usecprj) 283 type(paw_ij_type),intent(in) :: paw_ij(my_natom*psps%usepaw),paw_ij1(my_natom*psps%usepaw) 284 type(pawfgrtab_type),intent(inout) :: pawfgrtab(my_natom*psps%usepaw) 285 type(pawrhoij_type),intent(in) :: pawrhoij(my_natom*psps%usepaw) 286 type(pawrhoij_type),target,intent(inout) :: pawrhoij1(my_natom*psps%usepaw) 287 type(pawtab_type), intent(in) :: pawtab(ntypat*psps%usepaw) 288 type(wfk_t),intent(inout) :: ddk_f(4) 289 290 !Local variables------------------------------- 291 !scalars 292 integer,parameter :: level=13 293 integer :: bd2tot_index,bdtot_index,buffer_size,counter,cplex_rhoij 294 integer :: iband,nlines_done,ibdkpt,ibg,ibg1,ibgq,icg,icg1,icgq,ierr 295 integer :: ii,ikg,ikg1,ikpt,ilm,index1,ispden,iscf_mod,isppol,istwf_k 296 integer :: mbd2kpsp,mbdkpsp,mcgq,mcgq_disk,mcprjq 297 integer :: mcprjq_disk,me,n1,n2,n3,n4,n5,n6,nband_k,nband_kq,nkpg,nkpg1 298 integer :: nband_eff 299 integer :: nnsclo_now,npw1_k,npw_k,nspden_rhoij,qphase_rhoij,spaceworld,test_dot 300 integer :: nband_me 301 logical :: has_vectornd,has_vxctau,paral_atom,qne0 302 real(dp) :: arg,wtk_k 303 type(gs_hamiltonian_type) :: gs_hamkq 304 type(rf_hamiltonian_type) :: rf_hamkq,rf_hamk_dir2 305 !arrays 306 integer,allocatable :: kg1_k(:,:),kg_k(:,:) 307 integer, pointer :: my_atmtab(:) 308 real(dp) :: kpoint(3),kpq(3) 309 real(dp) :: tsec(2) 310 real(dp),allocatable :: buffer1(:) 311 real(dp),allocatable :: ddkinpw(:),dkinpw(:),dkinpw2(:) 312 real(dp),allocatable :: doccde_k(:),doccde_kq(:) 313 real(dp),allocatable :: edocc_k(:),eeig0_k(:),eig0_k(:),eig0_kq(:),eig1_k(:) 314 real(dp),allocatable :: ek0_k(:),ek1_k(:),eloc0_k(:),end0_k(:),end1_k(:),enl0_k(:),enl1_k(:) 315 real(dp),allocatable :: evxctau0_k(:),evxctau1_k(:) 316 real(dp),allocatable :: ffnl1(:,:,:,:),ffnl1_test(:,:,:,:),ffnlk(:,:,:,:) 317 real(dp),allocatable :: grad_berry(:,:,:),kinpw1(:),kpg1_k(:,:) 318 real(dp),allocatable :: kpg_k(:,:),occ_k(:),occ_kq(:) 319 real(dp),allocatable :: ph3d(:,:,:),ph3d1(:,:,:),resid_k(:) 320 real(dp),allocatable :: rho1wfg(:,:),rho1wfr(:,:),rhoaug1(:,:,:,:),rocceig(:,:) 321 real(dp),allocatable :: vectornd_pac(:,:,:,:,:),vectornd_pac_idir(:,:,:,:),vlocal(:,:,:,:),vlocal1(:,:,:,:) 322 real(dp),allocatable :: vxctaulocal(:,:,:,:,:) 323 real(dp),allocatable :: ylm1_k(:,:),ylm_k(:,:),ylmgr1_k(:,:,:) 324 type(pawrhoij_type),pointer :: pawrhoij1_unsym(:) 325 326 ! ********************************************************************* 327 328 DBG_ENTER('COLL') 329 330 !Keep track of total time spent in this routine 331 call timab(121,1,tsec) 332 call timab(124,1,tsec) 333 334 !Retrieve parallelism data 335 spaceworld=mpi_enreg%comm_cell 336 me=mpi_enreg%me_kpt 337 paral_atom=(my_natom/=natom) 338 my_atmtab=>mpi_enreg%my_atmtab 339 340 if (dtset%berryopt== 4.or.dtset%berryopt== 6.or.dtset%berryopt== 7.or.& 341 & dtset%berryopt==14.or.dtset%berryopt==16.or.dtset%berryopt==17) then 342 ABI_MALLOC(grad_berry,(2,mpw1,dtefield%mband_occ)) 343 else 344 ABI_MALLOC(grad_berry,(0,0,0)) 345 end if 346 347 !Test size of FFT grids (1 grid in norm-conserving, 2 grids in PAW) 348 if ((psps%usepaw==1.and.pawfgr%nfft/=nfftf).or.(psps%usepaw==0.and.dtset%nfft/=nfftf)) then 349 ABI_BUG('wrong values for nfft, nfftf!') 350 end if 351 352 !The value of iscf must be modified if ddk perturbation, see dfpt_looppert.f 353 iscf_mod=dtset%iscf;if(ipert==natom+1.or.ipert==natom+10.or.ipert==natom+11) iscf_mod=-3 354 355 edocc=zero ; eeig0=zero ; ek0=zero ; ek1=zero 356 eloc0=zero ; end0=zero ; end1=zero ; enl0=zero ; enl1=zero 357 evxctau0=zero; evxctau1=zero 358 bdtot_index=0 359 bd2tot_index=0 360 ibg=0;icg=0 361 ibgq=0;icgq=0 362 ibg1=0;icg1=0 363 mbdkpsp=mband*nkpt_rbz*nsppol 364 mbd2kpsp=2*mband**2*nkpt_rbz*nsppol 365 366 n1=dtset%ngfft(1); n2=dtset%ngfft(2); n3=dtset%ngfft(3) 367 n4=dtset%ngfft(4); n5=dtset%ngfft(5); n6=dtset%ngfft(6) 368 qne0=(qphon(1)**2+qphon(2)**2+qphon(3)**2>=tol14) 369 370 !Initialize PW 1st-order density if needed 371 !Also store old rho1 in case of density mixing 372 if (iscf_mod>0) then 373 if (optres==1) nvresid1=rhor1 374 if (psps%usepaw==0) then 375 rhor1(:,:)=zero 376 else 377 ABI_MALLOC(rho1wfr,(cplex*dtset%nfft,dtset%nspden)) 378 ABI_MALLOC(rho1wfg,(2,dtset%nfft)) 379 rho1wfr(:,:)=zero 380 end if 381 end if 382 383 !Set max number of non-self-consistent loops nnsclo_now for use in dfpt_vtowfk 384 if(iscf_mod<=0 .and. iscf_mod/=-3)then 385 nnsclo_now=dtset%nstep 386 else 387 if(dtset%nnsclo>0)then 388 nnsclo_now=dtset%nnsclo 389 else 390 nnsclo_now=1 391 end if 392 if(dbl_nnsclo==1) nnsclo_now=nnsclo_now*2 393 end if 394 395 !Prepare GS k+q wf 396 mcgq=mpw1*dtset%nspinor*mband_mem*mkqmem*nsppol;mcgq_disk=0 397 398 !Prepare RF PAW files 399 if (psps%usepaw==1) then 400 mcprjq=dtset%nspinor*mband_mem*mkqmem*nsppol*usecprj;mcprjq_disk=0 401 else 402 mcprjq=0;mcprjq_disk=0 403 end if 404 405 !Initialisation of the wfdot file in case of electric field (or 2nd order Sternheimer equation) 406 test_dot=0 407 if (ipert==natom+2.and.sum((dtset%qptn(1:3))**2 )<=tol7.and.& 408 & (dtset%berryopt/= 4.and.dtset%berryopt/= 6.and.dtset%berryopt/= 7.and.& 409 & dtset%berryopt/=14.and.dtset%berryopt/=16.and.dtset%berryopt/=17).or.& 410 & (ipert==natom+10.or.ipert==natom+11)) then 411 test_dot=1 412 end if 413 414 !==== Initialize most of the Hamiltonian (and derivative) ==== 415 !1) Allocate all arrays and initialize quantities that do not depend on k and spin. 416 !2) Perform the setup needed for the non-local factors: 417 !* Norm-conserving: Constant kleimann-Bylander energies are copied from psps to gs_hamk. 418 !* PAW: Initialize the overlap coefficients and allocate the Dij coefficients. 419 420 call init_hamiltonian(gs_hamkq,psps,pawtab,dtset%nspinor,nsppol,nspden,natom,& 421 & dtset%typat,xred,dtset%nfft,dtset%mgfft,dtset%ngfft,rprimd,dtset%nloalg,& 422 & paw_ij=paw_ij,comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab,mpi_spintab=mpi_enreg%my_isppoltab,& 423 & usecprj=usecprj,ph1d=ph1d,nucdipmom=dtset%nucdipmom,gpu_option=dtset%gpu_option) 424 425 call init_rf_hamiltonian(cplex,gs_hamkq,ipert,rf_hamkq,has_e1kbsc=.true.,paw_ij1=paw_ij1,& 426 & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab,mpi_spintab=mpi_enreg%my_isppoltab) 427 if ((ipert==natom+10.and.idir>3).or.ipert==natom+11) then 428 call init_rf_hamiltonian(cplex,gs_hamkq,ipert,rf_hamk_dir2,has_e1kbsc=.true.,paw_ij1=paw_ij1,& 429 & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab,mpi_spintab=mpi_enreg%my_isppoltab) 430 end if 431 432 !PAW:allocate memory for non-symetrized 1st-order occupancies matrix (pawrhoij1) 433 pawrhoij1_unsym => pawrhoij1 434 if (psps%usepaw==1.and.iscf_mod>0) then 435 if (paral_atom) then 436 ABI_MALLOC(pawrhoij1_unsym,(natom)) 437 call pawrhoij_inquire_dim(cplex_rhoij=cplex_rhoij,qphase_rhoij=qphase_rhoij,nspden_rhoij=nspden_rhoij,& 438 & nspden=dtset%nspden,spnorb=dtset%pawspnorb,cplex=cplex,cpxocc=dtset%pawcpxocc) 439 call pawrhoij_alloc(pawrhoij1_unsym,cplex_rhoij,nspden_rhoij,dtset%nspinor,& 440 & dtset%nsppol,dtset%typat,qphase=qphase_rhoij,pawtab=pawtab,& 441 & use_rhoijp=0,use_rhoij_=1) 442 else 443 pawrhoij1_unsym => pawrhoij1 444 call pawrhoij_init_unpacked(pawrhoij1_unsym) 445 end if 446 end if 447 448 ABI_MALLOC(rhoaug1,(cplex*n4,n5,n6,gs_hamkq%nvloc)) 449 ABI_MALLOC(vlocal,(n4,n5,n6,gs_hamkq%nvloc)) 450 ABI_MALLOC(vlocal1,(cplex*n4,n5,n6,gs_hamkq%nvloc)) 451 452 has_vxctau = ( size(vxctau) > 0 ) 453 if(has_vxctau) then 454 ABI_MALLOC(vxctaulocal,(n4,n5,n6,gs_hamkq%nvloc,4)) 455 end if 456 457 has_vectornd = (with_vectornd .EQ. 1) 458 if(has_vectornd) then 459 ABI_MALLOC(vectornd_pac,(n4,n5,n6,gs_hamkq%nvloc,3)) 460 ABI_MALLOC(vectornd_pac_idir,(n4,n5,n6,gs_hamkq%nvloc)) 461 end if 462 463 nlines_done = 0 464 resid = zero 465 466 !LOOP OVER SPINS 467 do isppol=1,nsppol 468 469 ! Rewind kpgsph data file if needed: 470 ikg=0;ikg1=0 471 472 ! Set up local potential vlocal1 with proper dimensioning, from vtrial1 473 ! Same thing for vlocal from vtrial Also take into account the spin. 474 475 call rf_transgrid_and_pack(isppol,nspden,psps%usepaw,cplex,nfftf,dtset%nfft,dtset%ngfft,& 476 & gs_hamkq%nvloc,pawfgr,mpi_enreg,vtrial,vtrial1,vlocal,vlocal1) 477 478 ! Continue to initialize the Hamiltonian 479 call gs_hamkq%load_spin(isppol,vlocal=vlocal,with_nonlocal=.true.) 480 call rf_hamkq%load_spin(isppol,vlocal1=vlocal1,with_nonlocal=.true.) 481 if ((ipert==natom+10.and.idir>3).or.ipert==natom+11) then 482 call rf_hamk_dir2%load_spin(isppol,with_nonlocal=.true.) 483 if (ipert==natom+11) then ! load vlocal1 484 call rf_hamk_dir2%load_spin(isppol,vlocal1=vlocal1) 485 end if 486 end if 487 488 if (ipert==natom+5) then !SPr deb, in case of magnetic field perturbation, no non-local 489 call rf_hamkq%load_spin(isppol,vlocal1=vlocal1) 490 end if 491 492 ! Nullify contribution to 1st-order density from this k-point 493 rhoaug1(:,:,:,:)=zero 494 495 ! if vectornd is present, set it up for addition to gs_hamkq and rf_hamkq. 496 ! Note that it must be done for the three Cartesian directions. Also, the following 497 ! code assumes explicitly and implicitly that nvloc = 1. This should eventually be generalized. 498 if(has_vectornd) then 499 call gspot_transgrid_and_pack(isppol, psps%usepaw, dtset%paral_kgb, dtset%nfft, dtset%ngfft, nfftf, & 500 & dtset%nspden, gs_hamkq%nvloc, 3, pawfgr, mpi_enreg, vectornd, vectornd_pac) 501 call gs_hamkq%load_spin(isppol, vectornd=vectornd_pac) 502 vectornd_pac_idir(:,:,:,:)=vectornd_pac(:,:,:,:,idir) 503 call rf_hamkq%load_spin(isppol, vectornd=vectornd_pac_idir) 504 end if 505 506 !! add vxctau for mGGA to GS hamiltonian and RF hamiltonian 507 if (has_vxctau) then 508 call gspot_transgrid_and_pack(isppol, psps%usepaw, dtset%paral_kgb, dtset%nfft, dtset%ngfft, nfftf, & 509 dtset%nspden, gs_hamkq%nvloc, 4, pawfgr, mpi_enreg, vxctau, vxctaulocal) 510 call gs_hamkq%load_spin(isppol, vxctaulocal=vxctaulocal) 511 call rf_hamkq%load_spin(isppol, vxctaulocal=vxctaulocal) 512 end if 513 514 call timab(125,1,tsec) 515 516 !====================================================================== 517 !============== BIG FAT K POINT LOOP ================================ 518 !====================================================================== 519 520 do ikpt=1,nkpt_rbz 521 counter=100*ikpt+isppol 522 523 nband_k = nband_rbz(ikpt+(isppol-1)*nkpt_rbz) 524 ! enables variable nband less than the block size in serial case 525 nband_me = proc_distrb_nband(mpi_enreg%proc_distrb,ikpt,nband_k,isppol,me) 526 istwf_k = istwfk_rbz(ikpt) 527 npw_k = npwarr(ikpt,1) 528 npw1_k = npwar1(ikpt,1) 529 wtk_k = wtk_rbz(ikpt) 530 531 if(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,isppol,me)) then 532 eigen1(1+bd2tot_index : 2*nband_k**2+bd2tot_index) = zero 533 resid(1+bdtot_index : nband_k+bdtot_index) = zero 534 bdtot_index=bdtot_index+nband_k 535 bd2tot_index=bd2tot_index+2*nband_k**2 536 537 cycle ! Skip the rest of the k-point loop 538 end if 539 540 kpoint(:)=kpt_rbz(:,ikpt) 541 kpq(:)=kpoint(:);if (ipert<natom+3.or.ipert==natom+5) kpq(:)=kpq(:)+qphon(1:3) 542 ABI_MALLOC(kg_k,(3,npw_k)) 543 ABI_MALLOC(kg1_k,(3,npw1_k)) 544 ABI_MALLOC(ylm_k,(npw_k,psps%mpsang*psps%mpsang*psps%useylm)) 545 ABI_MALLOC(ylm1_k,(npw1_k,psps%mpsang*psps%mpsang*psps%useylm)) 546 ABI_MALLOC(ylmgr1_k,(npw1_k,3+6*((ipert-natom)/10),psps%mpsang*psps%mpsang*psps%useylm*useylmgr1)) 547 ABI_MALLOC(doccde_k,(nband_k)) 548 ABI_MALLOC(doccde_kq,(nband_k)) 549 ABI_MALLOC(eig0_k,(nband_k)) 550 ABI_MALLOC(eig0_kq,(nband_k)) 551 ABI_MALLOC(eig1_k,(2*nband_k**2)) 552 ABI_MALLOC(edocc_k,(nband_k)) 553 ABI_MALLOC(eeig0_k,(nband_k)) 554 ABI_MALLOC(ek0_k,(nband_k)) 555 ABI_MALLOC(ek1_k,(nband_k)) 556 ABI_MALLOC(eloc0_k,(nband_k)) 557 ABI_MALLOC(end0_k,(nband_k)) 558 ABI_MALLOC(end1_k,(nband_k)) 559 ABI_MALLOC(enl0_k,(nband_k)) 560 ABI_MALLOC(enl1_k,(nband_k)) 561 ABI_MALLOC(evxctau0_k,(nband_k)) 562 ABI_MALLOC(evxctau1_k,(nband_k)) 563 ABI_MALLOC(occ_k,(nband_k)) 564 ABI_MALLOC(occ_kq,(nband_k)) 565 ABI_MALLOC(resid_k,(nband_k)) 566 ABI_MALLOC(rocceig,(nband_k,nband_k)) 567 568 eig1_k(:)=zero 569 eig0_k(:)=eigen0(1+bdtot_index:nband_k+bdtot_index) 570 eig0_kq(:)=eigenq(1+bdtot_index:nband_k+bdtot_index) 571 edocc_k(:)=zero 572 eeig0_k(:)=zero ; ek0_k(:)=zero ; ek1_k(:)=zero 573 eloc0_k(:)=zero ; end0_k(:)=zero ; end1_k(:)=zero 574 enl0_k(:)=zero ; enl1_k(:)=zero 575 evxctau0_k(:)=zero; evxctau1_k(:)=zero 576 occ_k(:)=occ_rbz(1+bdtot_index:nband_k+bdtot_index) 577 occ_kq(:)=occkq(1+bdtot_index:nband_k+bdtot_index) 578 doccde_k(:)=doccde_rbz(1+bdtot_index:nband_k+bdtot_index) 579 doccde_kq(:)=docckqde(1+bdtot_index:nband_k+bdtot_index) 580 resid_k(:)=zero 581 582 ! For each pair of active bands (m,n), generates the ratios 583 ! rocceig(m,n)=(occ_kq(m)-occ_k(n))/(eig0_kq(m)-eig0_k(n)) 584 ! and decide to which band to attribute it. 585 call occeig(doccde_k,doccde_kq,eig0_k,eig0_kq,nband_k,dtset%occopt,occ_k,occ_kq,rocceig) 586 587 ! These arrays are not needed anymore. 588 ABI_FREE(doccde_k) 589 ABI_FREE(doccde_kq) 590 ABI_FREE(occ_kq) 591 592 kg_k(:,1:npw_k)=kg(:,1+ikg:npw_k+ikg) 593 if (psps%useylm==1) then 594 do ilm=1,psps%mpsang*psps%mpsang 595 ylm_k(1:npw_k,ilm)=ylm(1+ikg:npw_k+ikg,ilm) 596 end do 597 end if 598 599 ! Get (k+q+G) wave vectors and associated spherical harmonics 600 kg1_k(:,1:npw1_k)=kg1(:,1+ikg1:npw1_k+ikg1) 601 if (psps%useylm==1) then 602 do ilm=1,psps%mpsang*psps%mpsang 603 ylm1_k(1:npw1_k,ilm)=ylm1(1+ikg1:npw1_k+ikg1,ilm) 604 end do 605 if (useylmgr1==1) then 606 do ilm=1,psps%mpsang*psps%mpsang 607 do ii=1,3+6*((ipert-natom)/10) 608 ylmgr1_k(1:npw1_k,ii,ilm)=ylmgr1(1+ikg1:npw1_k+ikg1,ii,ilm) 609 end do 610 end do 611 end if 612 end if 613 614 ! Set up the ground-state Hamiltonian, and some parts of the 1st-order Hamiltonian 615 call getgh1c_setup(gs_hamkq,rf_hamkq,dtset,psps,& ! In 616 kpoint,kpq,idir,ipert,natom,rmet,gprimd,gmet,istwf_k,& ! In 617 npw_k,npw1_k,useylmgr1,kg_k,ylm_k,kg1_k,ylm1_k,ylmgr1_k,& ! In 618 dkinpw,nkpg,nkpg1,kpg_k,kpg1_k,kinpw1,ffnlk,ffnl1,ph3d,ph3d1,& ! Out 619 ddkinpw=ddkinpw,dkinpw2=dkinpw2,rf_hamk_dir2=rf_hamk_dir2,& ! Optional 620 ffnl1_test=ffnl1_test) ! Optional 621 622 ! Compute the gradient of the Berry-phase term 623 if (dtset%berryopt== 4.or.dtset%berryopt== 6.or.dtset%berryopt== 7.or.& 624 & dtset%berryopt==14.or.dtset%berryopt==16.or.dtset%berryopt==17) then 625 if (ipert<=natom) then 626 ! phonon perturbation 627 call dfptff_gradberry(cg,cg1,dtefield,grad_berry,ikpt,isppol,mband,mband_mem,mpw,mpw1,mkmem,mk1mem,& 628 & mpi_enreg,nkpt_rbz,& 629 & npwarr,npwar1,dtset%nspinor,nsppol,qmat,pwindall) 630 else 631 ! electric field perturbation 632 call dfptff_gbefd(cg,cg1,dtefield,grad_berry,idir,ikpt,isppol,mband,mband_mem,mpw,mpw1,mkmem,mk1mem,& 633 & mpi_enreg,nkpt_rbz,& 634 & npwarr,npwar1,dtset%nspinor,& 635 & nsppol,qmat,pwindall,rprimd) 636 end if 637 end if 638 639 ! Free some memory before calling dfpt_vtowfk 640 ABI_FREE(ylm_k) 641 ABI_FREE(ylm1_k) 642 ABI_FREE(ylmgr1_k) 643 644 ! Compute the eigenvalues, wavefunction, residuals, 645 ! contributions to kinetic energy, nonlocal energy, forces, 646 ! and update of 1st-order density to this k-point and this spin polarization. 647 nband_kq = nband_k !Note that the calculation only works for same number of bands on all K points. 648 ! Note that dfpt_vtowfk is called with kpoint, while kpt is used inside vtowfk3 649 call dfpt_vtowfk(cg,cgq,cg1,cg1_active,cplex,cprj,cprjq,cprj1,dim_eig2rf,dtfil,& 650 & dtset,edocc_k,eeig0_k,eig0_k,eig0_kq,eig1_k,ek0_k,ek1_k,eloc0_k,end0_k,end1_k,enl0_k,enl1_k,& 651 & evxctau0_k,evxctau1_k,fermie1,& 652 & ffnl1,ffnl1_test,gh0c1_set,gh1c_set,grad_berry,gs_hamkq,ibg,ibgq,ibg1,icg,icgq,icg1,idir,ikpt,ipert,isppol,& 653 & mband,mband_mem,mcgq,mcprjq,mkmem,mk1mem,mpi_enreg,mpw,mpw1,natom,nband_k,ncpgr,nnsclo_now,& 654 & npw_k,npw1_k,dtset%nspinor,nsppol,n4,n5,n6,occ_k,pawrhoij1_unsym,prtvol,psps,resid_k,& 655 & rf_hamkq,rf_hamk_dir2,rhoaug1,rocceig,ddk_f,wtk_k,nlines_done,cg1_out) 656 657 ! Free temporary storage 658 ABI_FREE(kinpw1) 659 ABI_FREE(kg_k) 660 ABI_FREE(kg1_k) 661 ABI_FREE(kpg_k) 662 ABI_FREE(kpg1_k) 663 ABI_FREE(dkinpw) 664 if (ipert==natom+10) then 665 ABI_FREE(ddkinpw) 666 if (idir>3) then 667 ABI_FREE(dkinpw2) 668 end if 669 end if 670 ABI_FREE(ffnlk) 671 ABI_FREE(ffnl1) 672 if (allocated(ffnl1_test)) then 673 ABI_FREE(ffnl1_test) 674 end if 675 ABI_FREE(eig0_k) 676 ABI_FREE(eig0_kq) 677 ABI_FREE(rocceig) 678 ABI_FREE(ph3d) 679 if (allocated(ph3d1)) then 680 ABI_FREE(ph3d1) 681 end if 682 683 ! Save eigenvalues (hartree), residuals (hartree**2) 684 eigen1 (1+bd2tot_index : 2*nband_k**2+bd2tot_index) = eig1_k(:) 685 resid (1+bdtot_index : nband_k+bdtot_index) = resid_k(:) 686 687 ! Accumulate sum over k points for nonlocal and kinetic energies, 688 ! also accumulate gradients of Enonlocal: 689 if (iscf_mod>0 .or. iscf_mod==-3 .or. iscf_mod==-2)then 690 do iband=1,nband_k 691 edocc=edocc+wtk_k*occ_k(iband)*edocc_k(iband) 692 eeig0=eeig0+wtk_k*occ_k(iband)*eeig0_k(iband) 693 ek0=ek0+wtk_k*occ_k(iband)*ek0_k(iband) 694 ek1=ek1+wtk_k*occ_k(iband)*ek1_k(iband) 695 eloc0=eloc0+wtk_k*occ_k(iband)*eloc0_k(iband) 696 end0=end0+wtk_k*occ_k(iband)*end0_k(iband) 697 end1=end1+wtk_k*occ_k(iband)*end1_k(iband) 698 enl0=enl0+wtk_k*occ_k(iband)*enl0_k(iband) 699 enl1=enl1+wtk_k*occ_k(iband)*enl1_k(iband) 700 evxctau0=evxctau0+wtk_k*occ_k(iband)*evxctau0_k(iband) 701 evxctau1=evxctau1+wtk_k*occ_k(iband)*evxctau1_k(iband) 702 end do 703 end if 704 705 ABI_FREE(eig1_k) 706 ABI_FREE(occ_k) 707 ABI_FREE(resid_k) 708 ABI_FREE(edocc_k) 709 ABI_FREE(eeig0_k) 710 ABI_FREE(ek0_k) 711 ABI_FREE(ek1_k) 712 ABI_FREE(eloc0_k) 713 ABI_FREE(end0_k) 714 ABI_FREE(end1_k) 715 ABI_FREE(enl0_k) 716 ABI_FREE(enl1_k) 717 ABI_FREE(evxctau0_k) 718 ABI_FREE(evxctau1_k) 719 720 ! Keep track of total number of bands (all k points so far, even for k points not treated by me) 721 bdtot_index=bdtot_index+nband_k 722 bd2tot_index=bd2tot_index+2*nband_k**2 723 724 ! Shift array memory 725 if (mkmem/=0) then 726 ibg=ibg+dtset%nspinor*nband_me 727 icg=icg+npw_k*dtset%nspinor*nband_me 728 ikg=ikg+npw_k 729 end if 730 if (mkqmem/=0) then 731 ibgq=ibgq+dtset%nspinor*nband_me 732 icgq=icgq+npw1_k*dtset%nspinor*nband_me 733 end if 734 if (mk1mem/=0) then 735 ibg1=ibg1+dtset%nspinor*nband_me 736 icg1=icg1+npw1_k*dtset%nspinor*nband_me 737 ikg1=ikg1+npw1_k 738 end if 739 740 end do !ikpt loop 741 742 !====================================================================== 743 !================== END BIG K POINT LOOP ============================ 744 !====================================================================== 745 746 call timab(125,2,tsec) 747 748 ! Transfer density on augmented fft grid to normal fft grid in real space. Also take into account the spin. 749 ! FR EB for the non-collinear part see vtorho.F90 750 if(iscf_mod>0) then 751 if (psps%usepaw==0) then 752 call fftpac(isppol,mpi_enreg,nspden,cplex*n1,n2,n3,cplex*n4,n5,n6,dtset%ngfft,rhor1,rhoaug1(:,:,:,1),1) 753 if(nspden==4)then 754 do ispden=2,4 755 call fftpac(ispden,mpi_enreg,nspden,cplex*n1,n2,n3,cplex*n4,n5,n6,dtset%ngfft,rhor1,rhoaug1(:,:,:,ispden),1) 756 end do 757 end if 758 else 759 call fftpac(isppol,mpi_enreg,nspden,cplex*n1,n2,n3,cplex*n4,n5,n6,dtset%ngfft,rho1wfr,rhoaug1(:,:,:,1),1) 760 if(nspden==4)then 761 do ispden=2,4 762 call fftpac(ispden,mpi_enreg,nspden,cplex*n1,n2,n3,cplex*n4,n5,n6,dtset%ngfft,rho1wfr,rhoaug1(:,:,:,ispden),1) 763 end do 764 end if 765 end if 766 end if 767 768 end do ! End loop over spins 769 770 !More memory cleaning 771 call gs_hamkq%free() 772 call rf_hamkq%free() 773 if ((ipert==natom+10.and.idir>3).or.ipert==natom+11) then 774 call rf_hamk_dir2%free() 775 end if 776 ABI_FREE(rhoaug1) 777 ABI_FREE(vlocal) 778 ABI_FREE(vlocal1) 779 if(allocated(vxctaulocal)) then 780 ABI_FREE(vxctaulocal) 781 end if 782 if(allocated(vectornd_pac)) then 783 ABI_FREE(vectornd_pac) 784 end if 785 if(allocated(vectornd_pac_idir)) then 786 ABI_FREE(vectornd_pac_idir) 787 end if 788 789 call timab(124,2,tsec) 790 791 !=== MPI communications ================== 792 if(xmpi_paral==1)then 793 call timab(129,1,tsec) 794 795 ! MG: For the record, buffer1 can be pretty big if we have dense k-meshes e.g. metals 796 ! this is what you get with nband 26 and ngkpt 42**3: 797 ! 798 ! [0] <var=buffer1, A@m_dfpt_vtorho.F90:773, addr=0x150d5a59e010, size_mb=401.130> 799 ! 800 ! and this can lead to OOM if we have 2Gb per core also because xmpi_sum allocates another array of the same size! 801 ! TODO: Avoid packing rhor1 in buffer 802 803 ! Compute buffer size 804 buffer_size=11+mbd2kpsp+mbdkpsp 805 if (iscf_mod>0) then 806 buffer_size=buffer_size+cplex*dtset%nfft*nspden 807 end if 808 ABI_MALLOC(buffer1,(buffer_size)) 809 810 ! Pack rhor1,edocc,eeig0,ek0,ek1,eloc0,end0,end1,enl0,enl1,evxctau0,evxctau1,eigen1,resid 811 if (iscf_mod>0) then 812 index1=cplex*dtset%nfft*nspden 813 if (psps%usepaw==0) then 814 buffer1(1:index1)=reshape(rhor1 ,(/index1/)) 815 else 816 buffer1(1:index1)=reshape(rho1wfr,(/index1/)) 817 end if 818 else 819 index1=0 820 end if 821 buffer1(index1+1)=edocc;buffer1(index1+2)=eeig0 822 buffer1(index1+3)=ek0 ;buffer1(index1+4)=ek1 823 buffer1(index1+5)=eloc0;buffer1(index1+6)=enl0 824 buffer1(index1+7)=enl1 825 buffer1(index1+8)=end0;buffer1(index1+9)=end1 826 buffer1(index1+10)=evxctau0;buffer1(index1+11)=evxctau1 827 index1=index1+11 828 bdtot_index=0;bd2tot_index=0 829 do isppol=1,nsppol 830 do ikpt=1,nkpt_rbz 831 nband_k=nband_rbz(ikpt+(isppol-1)*nkpt_rbz) 832 buffer1(index1+1:index1+2*nband_k**2) = eigen1(bd2tot_index+1:bd2tot_index+2*nband_k**2) 833 buffer1(index1+2*nband_k**2+1:index1+2*nband_k**2+nband_k)= resid(bdtot_index+1:bdtot_index+nband_k) 834 bdtot_index=bdtot_index+nband_k 835 bd2tot_index=bd2tot_index+2*nband_k**2 836 index1=index1+2*nband_k**2+nband_k 837 end do 838 end do 839 if(index1<buffer_size)buffer1(index1+1:buffer_size)=zero 840 841 ! Build sum of everything 842 call timab(48,1,tsec) 843 call xmpi_sum(buffer1,buffer_size,spaceworld,ierr) 844 call timab(48,2,tsec) 845 846 ! Unpack the final result 847 if(iscf_mod>0) then 848 index1=cplex*dtset%nfft*nspden 849 if (psps%usepaw==0) then 850 rhor1(:,:) =reshape(buffer1(1:index1),(/cplex*dtset%nfft,nspden/)) 851 else 852 rho1wfr(:,:)=reshape(buffer1(1:index1),(/cplex*dtset%nfft,nspden/)) 853 end if 854 else 855 index1=0 856 end if 857 858 edocc=buffer1(index1+1);eeig0=buffer1(index1+2) 859 ek0=buffer1(index1+3) ;ek1=buffer1(index1+4) 860 eloc0=buffer1(index1+5);enl0=buffer1(index1+6) 861 enl1=buffer1(index1+7) 862 end0=buffer1(index1+8);end1=buffer1(index1+9) 863 evxctau0=buffer1(index1+10);evxctau1=buffer1(index1+11) 864 index1=index1+11 865 bdtot_index=0;bd2tot_index=0 866 do isppol=1,nsppol 867 do ikpt=1,nkpt_rbz 868 nband_k=nband_rbz(ikpt+(isppol-1)*nkpt_rbz) 869 eigen1(bd2tot_index+1:bd2tot_index+2*nband_k**2) = buffer1(index1+1:index1+2*nband_k**2) 870 resid(bdtot_index+1:bdtot_index+nband_k)= buffer1(index1+2*nband_k**2+1:index1+2*nband_k**2+nband_k) 871 bdtot_index=bdtot_index+nband_k 872 bd2tot_index=bd2tot_index+2*nband_k**2 873 index1=index1+2*nband_k**2+nband_k 874 end do 875 end do 876 ABI_FREE(buffer1) 877 878 ! Accumulate PAW occupancies 879 if (psps%usepaw==1.and.iscf_mod>0) then 880 call pawrhoij_mpisum_unpacked(pawrhoij1_unsym,spaceworld) 881 end if 882 883 call timab(129,2,tsec) 884 end if ! if kpt parallel 885 886 call timab(127,1,tsec) 887 888 !If needed, compute rhog1, and symmetrizes the density 889 if (iscf_mod > 0) then 890 891 ! In order to have the symrhg working in parallel on FFT coefficients, the size 892 ! of irzzon1 and phnons1 should be set to nfftot. Therefore, nsym\=1 does not work. 893 894 if(nspden==4) then 895 ! FR symrhg will manage correctly this rearrangement 896 rhor1(:,2)=rhor1(:,2)+(rhor1(:,1)+rhor1(:,4)) !(n+mx) 897 rhor1(:,3)=rhor1(:,3)+(rhor1(:,1)+rhor1(:,4)) !(n+my) 898 end if 899 ! 900 if (psps%usepaw==0) then 901 call symrhg(cplex,gprimd,irrzon1,mpi_enreg,dtset%nfft,dtset%nfft,dtset%ngfft,& 902 & nspden,nsppol,nsym1,phnons1,rhog1,rhor1,rprimd,symaf1,symrl1,tnons1) 903 else 904 call symrhg(cplex,gprimd,irrzon1,mpi_enreg,dtset%nfft,dtset%nfft,dtset%ngfft,& 905 & nspden,nsppol,nsym1,phnons1,rho1wfg,rho1wfr,rprimd,symaf1,symrl1,tnons1) 906 end if 907 ! We now have both rho(r) and rho(G), symmetrized, and if nsppol=2 908 ! we also have the spin-up density, symmetrized, in rhor1(:,2). 909 end if 910 911 ABI_FREE(grad_berry) 912 913 !Find largest residual over bands, k points, and spins except for nbdbuf highest bands 914 ibdkpt=1 915 residm=zero 916 do isppol=1,nsppol 917 do ikpt=1,nkpt_rbz 918 nband_k=nband_rbz(ikpt+(isppol-1)*nkpt_rbz) 919 nband_eff=max(1,nband_k-dtset%nbdbuf) 920 residm=max(residm,maxval(resid(ibdkpt:ibdkpt+nband_eff-1))) 921 ibdkpt=ibdkpt+nband_k 922 end do 923 end do 924 925 call timab(127,2,tsec) 926 927 if (iscf_mod>0) then 928 929 ! PAW: Build new 1st-order rhoij quantities then symetrize them 930 ! Compute and add the 1st-order compensation density to rho1wfr 931 ! to get the total 1st-order density 932 if (psps%usepaw==1) then 933 call pawmkrho(1,arg,cplex,gprimd,idir,indsy1,ipert,mpi_enreg,& 934 & my_natom,natom,nspden,nsym1,ntypat,dtset%paral_kgb,pawang,pawfgr,pawfgrtab,& 935 & dtset%pawprtvol,pawrhoij1,pawrhoij1_unsym,pawtab,dtset%qptn,rho1wfg,rho1wfr,& 936 & rhor1,rprimd,symaf1,symrc1,dtset%typat,ucvol,dtset%usewvl,xred,& 937 & pawang_sym=pawang1,pawnhat=nhat1,pawrhoij0=pawrhoij,rhog=rhog1) 938 ABI_FREE(rho1wfr) 939 ABI_FREE(rho1wfg) 940 if (paral_atom) then 941 call pawrhoij_free(pawrhoij1_unsym) 942 ABI_FREE(pawrhoij1_unsym) 943 end if 944 end if 945 946 ! Compute density residual (if required) and its squared norm 947 if (optres==1) then 948 nvresid1=rhor1-nvresid1 949 call sqnorm_v(1,nfftf,nres2,dtset%nspden,optres,nvresid1) 950 end if 951 end if ! iscf>0 952 953 call timab(121,2,tsec) 954 955 DBG_EXIT('COLL') 956 957 end subroutine dfpt_vtorho
ABINIT/m_dfpt_vtorho [ Modules ]
NAME
m_dfpt_vtorho
FUNCTION
COPYRIGHT
Copyright (C) 1998-2024 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 .
SOURCE
15 #if defined HAVE_CONFIG_H 16 #include "config.h" 17 #endif 18 19 #include "abi_common.h" 20 21 module m_dfpt_vtorho 22 23 use defs_basis 24 use m_abicore 25 use m_xmpi 26 use m_errors 27 use m_efield 28 use m_hamiltonian 29 use m_wfk 30 use m_cgtools 31 use m_dtset 32 use m_dtfil 33 34 35 use defs_datatypes, only : pseudopotential_type 36 use defs_abitypes, only : MPI_type 37 use m_time, only : timab 38 use m_occ, only : occeig 39 use m_hdr, only : hdr_skip, hdr_io 40 use m_pawang, only : pawang_type 41 use m_pawtab, only : pawtab_type 42 use m_paw_ij, only : paw_ij_type 43 use m_pawfgrtab,only : pawfgrtab_type 44 use m_pawrhoij, only : pawrhoij_type, pawrhoij_alloc, pawrhoij_free, & 45 & pawrhoij_init_unpacked, pawrhoij_free_unpacked, & 46 & pawrhoij_mpisum_unpacked, pawrhoij_inquire_dim 47 use m_pawcprj, only : pawcprj_type 48 use m_pawfgr, only : pawfgr_type 49 use m_paw_mkrho,only : pawmkrho 50 use m_fft, only : fftpac 51 use m_spacepar, only : symrhg 52 use m_getgh1c, only : rf_transgrid_and_pack, getgh1c_setup 53 use m_dfpt_vtowfk, only : dfpt_vtowfk 54 use m_dfpt_fef, only : dfptff_gradberry, dfptff_gbefd 55 use m_mpinfo, only : proc_distrb_cycle,proc_distrb_nband 56 use m_fourier_interpol, only : transgrid 57 58 implicit none 59 60 private