TABLE OF CONTENTS
ABINIT/dfpt_nselt [ Functions ]
NAME
dfpt_nselt
FUNCTION
This routine compute the non-stationary expression for the second derivative of the total energy, wrt strain for a whole row of mixed strain derivatives.
COPYRIGHT
Copyright (C) 1998-2018 ABINIT group (DRH, XG, DCA, GMR, MM, AR, MV, MB) 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 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 ecut=cut-off energy for plane wave basis sphere (Ha) ecutsm=smearing energy for plane wave kinetic energy (Ha) effmass_free=effective mass for electrons (1. in common case) gmet(3,3)=reciprocal space metric tensor in bohr**-2. gprimd(3,3)=dimensional reciprocal space primitive translations gsqcut=cutoff on (k+G)^2 (bohr^-2) idir=direction of the perturbation ipert=type of the perturbation 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 in the reduced BZ kxc(nfft,nkxc)=exchange and correlation kernel mband=maximum number of bands mgfft=maximum size of 1D FFTs mkmem =number of k points treated by this node. mk1mem =number of k points treated by this node (RF data). mpert =maximum number of ipert mpi_enreg=information about MPI parallelization mpsang= 1+maximum angular momentum for nonlocal pseudopotentials 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). maximum dimension for q points in grids for nonlocal form factors natom=number of atoms in cell. nband_rbz(nkpt_rbz*nsppol)=number of bands at each RF k point for each spin nfft=(effective) number of FFT grid points (for this processor) ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft nkpt_rbz=number of k points in the reduced BZ for this perturbation nkxc=second dimension of the kxc array. If /=0, the exchange-correlation kernel must be computed. nloalg(3)=governs the choice of the algorithm for non-local operator. 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 nspinor=number of spinorial components of the wavefunctions 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. occ_rbz(mband*nkpt_rbz*nsppol)=occupation number for each band and k in the reduced Brillouin zone (usually =2) ph1d(2,3*(2*mgfft+1)*natom)=one-dimensional structure factor information prtbbb=if 1, band-by-band decomposition (also dim of d2bbb) psps <type(pseudopotential_type)>=variables related to pseudopotentials qphon(3)=reduced coordinates for the phonon wavelength rhog(2,nfft)=array for Fourier transform of GS electron density rhor(nfft,nspden)=GS electron density in electrons/bohr**3. rhor1(cplex*nfft,nspden)=RF electron density in electrons/bohr**3. rmet(3,3)=real space metric (bohr**2) rprimd(3,3)=dimensional primitive translations in real space (bohr) symrc1(3,3,nsym1)=symmetry operations in reciprocal space typat(natom)=type integer for each atom in cell ucvol=unit cell volume in bohr**3. wtk_rbz(nkpt_rbz)=weight assigned to each k point in the reduced BZ xred(3,natom)=reduced dimensionless atomic coordinates ylm(mpw*mkmem,mpsang*mpsang)= real spherical harmonics for each G and k point ylm1(mpw1*mk1mem,mpsang*mpsang)= real spherical harmonics for each G and k+q point ylmgr(mpw*mkmem,3,mpsang*mpsang*useylm)= gradients of real spherical for each G and k point ylmgr1(mpw1*mk1mem,3,mpsang*mpsang*useylm)= gradients of real spherical for each G and k+g point
OUTPUT
blkflg(3,mpert,3,mpert)=flags for each element of the 2DTE (=1 if computed) d2bbb(2,3,3,mpert,mband,mband*prtbbb)=band by band decomposition of some second order derivatives d2lo(2,3,mpert,3,mpert)=local contributions to the 2DTEs d2nl(2,3,mpert,3,mpert)=non-local contributions to the 2DTEs
PARENTS
dfpt_scfcv
CHILDREN
destroy_hamiltonian,dfpt_mkcore,dfpt_mkvxcstr,dfpt_nsteltwf,dotprod_vn hartrestr,init_hamiltonian,stresssym,vlocalstr,wrtout,xmpi_barrier
SOURCE
101 #if defined HAVE_CONFIG_H 102 #include "config.h" 103 #endif 104 105 #include "abi_common.h" 106 107 108 subroutine dfpt_nselt(blkflg,cg,cg1,cplex,& 109 & d2bbb,d2lo,d2nl,ecut,ecutsm,effmass_free,& 110 & gmet,gprimd,gsqcut,idir,& 111 & ipert,istwfk_rbz,kg,kg1,kpt_rbz,kxc,mband,mgfft,& 112 & mkmem,mk1mem,mpert,mpi_enreg,mpsang,mpw,mpw1,& 113 & natom,nband_rbz,nfft,ngfft,& 114 & nkpt_rbz,nkxc,nloalg,npwarr,npwar1,nspden,nspinor,nsppol,& 115 & nsym1,ntypat,occ_rbz,& 116 & paral_kgb,ph1d,prtbbb,psps,qphon,rhog,& 117 & rhor,rhor1,rmet,rprimd,symrc1,typat,ucvol,& 118 & wtk_rbz,& 119 & xred,ylm,ylm1,ylmgr,ylmgr1) 120 121 use defs_basis 122 use defs_datatypes 123 use defs_abitypes 124 use m_errors 125 use m_profiling_abi 126 use m_xmpi 127 128 use m_cgtools, only : dotprod_vn 129 use m_pawtab, only : pawtab_type 130 use m_hamiltonian,only : init_hamiltonian,destroy_hamiltonian,gs_hamiltonian_type 131 132 !This section has been created automatically by the script Abilint (TD). 133 !Do not modify the following lines by hand. 134 #undef ABI_FUNC 135 #define ABI_FUNC 'dfpt_nselt' 136 use interfaces_14_hidewrite 137 use interfaces_32_util 138 use interfaces_41_geometry 139 use interfaces_72_response, except_this_one => dfpt_nselt 140 !End of the abilint section 141 142 implicit none 143 144 !Arguments ------------------------------- 145 !scalars 146 integer,intent(in) :: cplex,idir,ipert,mband,mgfft,mk1mem 147 integer,intent(in) :: mkmem,mpert,mpsang,mpw,mpw1,natom,nfft,nkpt_rbz 148 integer,intent(in) :: nkxc,nspden,nspinor,nsppol,nsym1,ntypat 149 integer,intent(in) :: paral_kgb,prtbbb 150 real(dp),intent(in) :: ecut,ecutsm,effmass_free,gsqcut,ucvol 151 type(MPI_type),intent(in) :: mpi_enreg 152 type(pseudopotential_type),intent(in) :: psps 153 !arrays 154 integer,intent(in) :: istwfk_rbz(nkpt_rbz) 155 integer,intent(in) :: kg(3,mpw*mkmem),kg1(3,mpw1*mk1mem) 156 integer,intent(in) :: nband_rbz(nkpt_rbz*nsppol),ngfft(18) 157 integer,intent(in) :: nloalg(3),npwar1(nkpt_rbz),npwarr(nkpt_rbz) 158 integer,intent(in) :: symrc1(3,3,nsym1),typat(natom) 159 integer,intent(inout) :: blkflg(3,mpert,3,mpert) 160 real(dp),intent(in) :: cg(2,mpw*nspinor*mband*mkmem*nsppol) 161 real(dp),intent(in) :: cg1(2,mpw1*nspinor*mband*mk1mem*nsppol) 162 real(dp),intent(in) :: gmet(3,3) 163 real(dp),intent(in) :: gprimd(3,3),kpt_rbz(3,nkpt_rbz),kxc(nfft,nkxc) 164 real(dp),intent(in) :: occ_rbz(mband*nkpt_rbz*nsppol) 165 real(dp),intent(in) :: ph1d(2,3*(2*mgfft+1)*natom),qphon(3),rhog(2,nfft) 166 real(dp),intent(in) :: rhor(nfft,nspden) 167 real(dp),intent(in) :: rhor1(cplex*nfft,nspden),rmet(3,3),rprimd(3,3) 168 real(dp),intent(in) :: wtk_rbz(nkpt_rbz),xred(3,natom) 169 real(dp),intent(in) :: ylm(mpw*mkmem,mpsang*mpsang*psps%useylm) 170 real(dp),intent(in) :: ylm1(mpw1*mk1mem,mpsang*mpsang*psps%useylm) 171 real(dp),intent(in) :: ylmgr(mpw*mkmem,3,mpsang*mpsang*psps%useylm) 172 real(dp),intent(in) :: ylmgr1(mpw1*mk1mem,3,mpsang*mpsang*psps%useylm) 173 real(dp),intent(out) :: d2bbb(2,3,3,mpert,mband,mband*prtbbb) 174 real(dp),intent(inout) :: d2lo(2,3,mpert,3,mpert) 175 real(dp),intent(inout) :: d2nl(2,3,mpert,3,mpert) 176 177 !Local variables------------------------------- 178 !scalars 179 integer :: ban2tot,bantot,bd2tot_index,bdtot_index 180 integer :: icg,icg1,idir1,ifft,ii,ikg,ikg1,ikpt,comm 181 integer :: ilm,ipert1,ispden,isppol,istr1,istwf_k 182 integer :: mbd2kpsp,mbdkpsp,me,n1,n2,n3,n3xccc,n4,n5,n6 183 integer :: nband_k,nfftot,npw1_k,npw_k,option 184 real(dp) :: doti,dotr 185 real(dp) :: wtk_k 186 character(len=500) :: message 187 type(gs_hamiltonian_type) :: gs_hamk 188 !arrays 189 integer :: ikpt_fbz(3) 190 integer,allocatable :: kg1_k(:,:),kg_k(:,:) 191 real(dp) :: kpoint(3),restr(6),dummy(0,0) 192 real(dp),allocatable :: d2bbb_k(:,:,:,:),d2nl_k(:,:,:) 193 real(dp),allocatable :: occ_k(:) 194 real(dp),allocatable :: vhartr01(:),vpsp1(:),vxc1(:,:),xccc3d1(:),ylm1_k(:,:) 195 real(dp),allocatable :: ylm_k(:,:),ylmgr1_k(:,:,:),ylmgr_k(:,:,:) 196 type(pawtab_type) :: pawtab_dum(0) 197 198 ! ********************************************************************* 199 200 !Init me 201 comm = mpi_enreg%comm_cell 202 me = mpi_enreg%me_kpt 203 204 !Unit numbers 205 206 !Zero only portion of nonlocal matrix to be computed here 207 d2nl(:,:,natom+3:natom+4,idir,ipert)=zero 208 bdtot_index=0 209 bd2tot_index=0 210 icg=0 211 icg1=0 212 mbdkpsp=mband*nkpt_rbz*nsppol 213 mbd2kpsp=2*mband**2*nkpt_rbz*nsppol 214 215 !Update list of computed matrix elements 216 if((ipert==natom+3) .or. (ipert==natom+4)) then 217 ! Eventually expand when strain coupling to other perturbations is 218 ! implemented 219 do ipert1=natom+3,natom+4 220 do idir1=1,3 221 blkflg(idir1,ipert1,idir,ipert)=1 222 end do 223 end do 224 end if 225 226 !allocate(enl1nk(mbdkpsp)) 227 ABI_ALLOCATE(d2bbb_k,(2,3,mband,mband*prtbbb)) 228 ABI_ALLOCATE(d2nl_k,(2,3,mpert)) 229 230 231 ABI_ALLOCATE(kg_k,(3,mpw)) 232 ABI_ALLOCATE(kg1_k,(3,mpw1)) 233 234 235 n1=ngfft(1) ; n2=ngfft(2) ; n3=ngfft(3) 236 n4=ngfft(4) ; n5=ngfft(5) ; n6=ngfft(6) 237 nfftot=n1*n2*n3 238 239 !Initialize Hamiltonian (k-independent terms) - NCPP only 240 call init_hamiltonian(gs_hamk,psps,pawtab_dum,nspinor,nsppol,nspden,natom,& 241 & typat,xred,nfft,mgfft,ngfft,rprimd,nloalg,ph1d=ph1d) 242 243 bantot = 0 244 ban2tot = 0 245 246 !LOOP OVER SPINS 247 do isppol=1,nsppol 248 249 if (nsppol/=1) then 250 write(message,*)' **** In dfpt_nselt for isppol=',isppol 251 call wrtout(std_out,message,'COLL') 252 end if 253 254 ikg=0 255 ikg1=0 256 257 ikpt_fbz(1:3)=0 258 259 ! BIG FAT k POINT LOOP 260 do ikpt=1,nkpt_rbz 261 262 nband_k=nband_rbz(ikpt+(isppol-1)*nkpt_rbz) 263 istwf_k=istwfk_rbz(ikpt) 264 npw_k=npwarr(ikpt) 265 npw1_k=npwar1(ikpt) 266 kpoint(:)=kpt_rbz(:,ikpt) 267 268 bantot = bantot + nband_k 269 ban2tot = ban2tot + 2*nband_k**2 270 271 if(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,isppol,me)) then 272 bdtot_index=bdtot_index+nband_k 273 bd2tot_index=bd2tot_index+2*nband_k**2 274 ! Skip the rest of the k-point loop 275 cycle 276 end if 277 278 ABI_ALLOCATE(occ_k,(nband_k)) 279 280 ABI_ALLOCATE(ylm_k,(npw_k,mpsang*mpsang)) 281 ABI_ALLOCATE(ylm1_k,(npw1_k,mpsang*mpsang)) 282 if (ipert==natom+3.or.ipert==natom+4) then 283 ABI_ALLOCATE(ylmgr_k,(npw_k,3,mpsang*mpsang)) 284 ABI_ALLOCATE(ylmgr1_k,(npw1_k,3,mpsang*mpsang)) 285 end if 286 287 ! enl1_k(:)=zero 288 d2nl_k(:,:,:)=zero 289 if(prtbbb==1)d2bbb_k(:,:,:,:)=zero 290 occ_k(:)=occ_rbz(1+bdtot_index:nband_k+bdtot_index) 291 292 kg_k(:,1:npw_k)=kg(:,1+ikg:npw_k+ikg) 293 if (psps%useylm==1) then 294 do ilm=1,mpsang*mpsang 295 ylm_k(1:npw_k,ilm)=ylm(1+ikg:npw_k+ikg,ilm) 296 end do 297 if (ipert==natom+3.or.ipert==natom+4) then 298 do ilm=1,mpsang*mpsang 299 do ii=1,3 300 ylmgr_k(1:npw_k,ii,ilm)=ylmgr(1+ikg:npw_k+ikg,ii,ilm) 301 end do 302 end do 303 end if 304 end if 305 306 wtk_k=wtk_rbz(ikpt) 307 308 kg1_k(:,:) = 0 309 310 kg1_k(:,1:npw1_k)=kg1(:,1+ikg1:npw1_k+ikg1) 311 if (psps%useylm==1) then 312 do ilm=1,mpsang*mpsang 313 ylm1_k(1:npw1_k,ilm)=ylm1(1+ikg1:npw1_k+ikg1,ilm) 314 end do 315 if (ipert==natom+3.or.ipert==natom+4) then 316 do ilm=1,mpsang*mpsang 317 do ii=1,3 318 ylmgr1_k(1:npw1_k,ii,ilm)=ylmgr1(1+ikg1:npw1_k+ikg1,ii,ilm) 319 end do 320 end do 321 end if 322 end if 323 324 ! Compute the eigenvalues, wavefunction, 325 ! contributions to kinetic energy, nonlocal energy, forces, 326 ! and update of rhor1 to this k-point and this spin polarization. 327 328 ! Note that dfpt_nsteltwf is called with kpoint, while kpt is used inside dfpt_vtowfk 329 call dfpt_nsteltwf(cg,cg1,d2nl_k,ecut,ecutsm,effmass_free,gs_hamk,icg,icg1,ikpt,isppol,& 330 & istwf_k,kg_k,kg1_k,kpoint,mband,mkmem,mk1mem,mpert,mpi_enreg,mpw,mpw1,natom,nband_k,& 331 & npw_k,npw1_k,nspinor,nsppol,ntypat,occ_k,psps,rmet,wtk_k,ylm_k,ylmgr_k) 332 d2nl(:,:,:,idir,ipert)=d2nl(:,:,:,idir,ipert)+d2nl_k(:,:,:) 333 if(prtbbb==1)then 334 d2bbb(:,:,idir,ipert,:,:) = d2bbb(:,:,idir,ipert,:,:) + & 335 & d2bbb_k(:,:,:,:) 336 end if 337 338 ABI_DEALLOCATE(occ_k) 339 340 ! Keep track of total number of bands (all k points so far, even for 341 ! k points not treated by me) 342 bdtot_index=bdtot_index+nband_k 343 bd2tot_index=bd2tot_index+2*nband_k**2 344 345 ! Shift array memory 346 if (mkmem/=0) then 347 icg=icg+npw_k*nspinor*nband_k 348 ikg=ikg+npw_k 349 end if 350 if (mk1mem/=0) then 351 icg1=icg1+npw1_k*nspinor*nband_k 352 ikg1=ikg1+npw1_k 353 end if 354 ABI_DEALLOCATE(ylm_k) 355 ABI_DEALLOCATE(ylm1_k) 356 if (ipert==natom+3.or.ipert==natom+4) then 357 ABI_DEALLOCATE(ylmgr_k) 358 ABI_DEALLOCATE(ylmgr1_k) 359 end if 360 361 end do ! End big k point loop 362 end do ! End loop over spins 363 364 if(xmpi_paral==1)then 365 call xmpi_barrier(comm) 366 call wrtout(std_out,' dfpt_nselt: loop on k-points and spins done in parallel','COLL') 367 end if 368 369 !Treat now varying occupation numbers 370 !if(occopt>=3 .and. occopt <=8) then 371 !SUPPRESSED metallic coding of vtorho 372 373 !Treat fixed occupation numbers 374 !else 375 376 !Accumulation over parallel processed now carried out for all terms 377 !in dfpt_nstdy.f 378 379 !End of test on varying or fixed occupation numbers 380 !end if 381 382 !The imaginary part of d2nl will be must be set to zero here since 383 !time-reversal symmetry will always be true for the strain peturbation. 384 !The symmetry-reduced kpt set will leave a non-zero imaginary part. 385 386 d2nl(2,:,natom+3:natom+4,idir,ipert)=zero 387 388 !Symmetrize the non-local contributions, 389 !as was needed for the stresses in a ground-state calculation 390 391 if (nsym1>1) then 392 ! Pack like symmetric-storage cartesian stress tensor 393 ii=0 394 do ipert1=natom+3,natom+4 395 do idir1=1,3 396 ii=ii+1 397 restr(ii)=d2nl(1,idir1,ipert1,idir,ipert) 398 end do 399 end do 400 ! Do the symmetrization using the ground state routine 401 call stresssym(gprimd,nsym1,restr,symrc1) 402 ! Unpack symmetrized stress tensor 403 ii=0 404 do ipert1=natom+3,natom+4 405 do idir1=1,3 406 ii=ii+1 407 d2nl(1,idir1,ipert1,idir,ipert)=restr(ii) 408 end do 409 end do 410 end if !nsym>1 411 412 !---------------------------------------------------------------------------- 413 !Now, treat the local contribution 414 415 ABI_ALLOCATE(vpsp1,(cplex*nfft)) 416 n3xccc=0 417 if(psps%n1xccc/=0)n3xccc=nfft 418 ABI_ALLOCATE(xccc3d1,(cplex*n3xccc)) 419 ABI_ALLOCATE(vxc1,(cplex*nfft,nspden)) 420 ABI_ALLOCATE(vhartr01,(nfft)) 421 xccc3d1(:)=zero 422 423 !Double loop over strain perturbations 424 do ipert1=natom+3,natom+4 425 do idir1=1,3 426 if(ipert1==natom+3) then 427 istr1=idir1 428 else 429 istr1=idir1+3 430 end if 431 432 ! Get first-order local potential. 433 call vlocalstr(gmet,gprimd,gsqcut,istr1,mgfft,mpi_enreg,& 434 & psps%mqgrid_vl,natom,gs_hamk%nattyp,nfft,ngfft,ntypat,paral_kgb,ph1d,psps%qgrid_vl,& 435 & ucvol,psps%vlspl,vpsp1) 436 437 ! Get first-order hartree potential. 438 call hartrestr(gsqcut,idir1,ipert1,mpi_enreg,natom,nfft,ngfft,& 439 & paral_kgb,rhog,rprimd,vhartr01) 440 441 ! Get first-order exchange-correlation potential 442 if(psps%n1xccc/=0)then 443 call dfpt_mkcore(cplex,idir1,ipert1,natom,ntypat,n1,psps%n1xccc,& 444 & n2,n3,qphon,rprimd,typat,ucvol,psps%xcccrc,psps%xccc1d,xccc3d1,xred) 445 end if ! psps%n1xccc/=0 446 447 option=0 448 call dfpt_mkvxcstr(cplex,idir1,ipert1,kxc,mpi_enreg,natom,nfft,ngfft,& 449 & dummy,dummy,nkxc,nspden,n3xccc,option,paral_kgb,qphon,rhor,rhor1,rprimd,& 450 & 0,0,vxc1,xccc3d1) 451 452 ! Combines density j2 with local potential j1 453 do ispden=1,min(nspden,2) 454 do ifft=1,cplex*nfft 455 vxc1(ifft,ispden)=vxc1(ifft,ispden)+vpsp1(ifft)+vhartr01(ifft) 456 end do 457 end do 458 call dotprod_vn(cplex,rhor1,dotr,doti,nfft,nfftot,nspden,2,vxc1,ucvol) 459 write(std_out,*) 460 d2lo(1,idir1,ipert1,idir,ipert)=dotr 461 d2lo(2,idir1,ipert1,idir,ipert)=doti 462 end do ! istr1 463 end do ! ipert1 464 465 call destroy_hamiltonian(gs_hamk) 466 467 ABI_DEALLOCATE(vxc1) 468 ABI_DEALLOCATE(xccc3d1) 469 ABI_DEALLOCATE(vhartr01) 470 471 ABI_DEALLOCATE(d2bbb_k) 472 ABI_DEALLOCATE(d2nl_k) 473 ABI_DEALLOCATE(kg_k) 474 ABI_DEALLOCATE(kg1_k) 475 ABI_DEALLOCATE(vpsp1) 476 477 end subroutine dfpt_nselt