TABLE OF CONTENTS
ABINIT/dfpt_nstwf [ Functions ]
NAME
dfpt_nstwf
FUNCTION
This routine computes the non-local contribution to the 2DTE matrix elements, in the non-stationary formulation Only for norm-conserving pseudopotentials (no PAW)
COPYRIGHT
Copyright (C) 1999-2018 ABINIT group (XG,AR,MB,MVer,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 at k
cg1(2,mpw1*nspinor*mband*mk1mem*nsppol)=pw coefficients of RF wavefunctions at k,q.
ddkfil(3)=unit numbers for the three possible ddk files for ipert1
equal to 0 if no dot file is available for this direction
dtset <type(dataset_type)>=all input variables for this dataset
eig_k(mband*nsppol)=GS eigenvalues at k (hartree)
eig1_k(2*nsppol*mband**2)=matrix of first-order eigenvalues (hartree)
gs_hamkq <type(gs_hamiltonian_type)>=all data for the Hamiltonian at k+q
icg=shift to be applied on the location of data in the array cg
icg1=shift to be applied on the location of data in the array cg1
idir=direction of the current perturbation
ikpt=number of the k-point
ipert=type of the perturbation
isppol=1 for unpolarized, 2 for spin-polarized
istwf_k=parameter that describes the storage of wfs
kg_k(3,npw_k)=reduced planewave coordinates.
kg1_k(3,npw1_k)=reduced planewave coordinates at k+q, with RF k points
kpt(3)=reduced coordinates of k point
kpq(3)=reduced coordinates of k+q point
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
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).
nband_k=number of bands at this k point for that spin polarization
npw_k=number of plane waves at this k point
npw1_k=number of plane waves at this k+q point
nsppol=1 for unpolarized, 2 for spin-polarized
occ_k(nband_k)=occupation number for each band (usually 2) for each k.
psps <type(pseudopotential_type)>=variables related to pseudopotentials
rmet(3,3)=real space metric (bohr**2)
ddks(3)<wfk_t>=struct info for for the three possible DDK files for ipert1
wtk_k=weight assigned to the k point.
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+q point
OUTPUT
d2bbb_k(2,3,mband,mband*prtbbb)=band by band decomposition of the second order derivatives, for the present k point, and perturbation idir, ipert d2nl_k(2,3,mpert)=non-local contributions to non-stationary 2DTE, for the present k point, and perturbation idir, ipert
TODO
XG 20141103 The localization tensor cannot be defined in the metallic case. It should not be computed.
PARENTS
dfpt_nstdy
CHILDREN
destroy_rf_hamiltonian,dotprod_g,gaugetransfo,getgh1c init_rf_hamiltonian,load_k_hamiltonian,load_k_rf_hamiltonian load_kprime_hamiltonian,mkffnl,mkkpg,timab,wfk_read_bks
SOURCE
75 #if defined HAVE_CONFIG_H 76 #include "config.h" 77 #endif 78 79 #include "abi_common.h" 80 81 subroutine dfpt_nstwf(cg,cg1,ddkfil,dtset,d2bbb_k,d2nl_k,eig_k,eig1_k,gs_hamkq,& 82 & icg,icg1,idir,ikpt,ipert,isppol,istwf_k,kg_k,kg1_k,kpt,kpq,& 83 & mkmem,mk1mem,mpert,mpi_enreg,mpw,mpw1,nband_k,npw_k,npw1_k,nsppol,& 84 & occ_k,psps,rmet,ddks,wtk_k,ylm,ylm1) 85 86 87 use defs_basis 88 use defs_datatypes 89 use defs_abitypes 90 use m_profiling_abi 91 use m_cgtools 92 use m_hamiltonian 93 use m_errors 94 use m_wfk 95 use m_xmpi 96 97 use m_pawcprj, only : pawcprj_type 98 99 !This section has been created automatically by the script Abilint (TD). 100 !Do not modify the following lines by hand. 101 #undef ABI_FUNC 102 #define ABI_FUNC 'dfpt_nstwf' 103 use interfaces_18_timing 104 use interfaces_66_nonlocal 105 use interfaces_66_wfs 106 use interfaces_72_response, except_this_one => dfpt_nstwf 107 !End of the abilint section 108 109 implicit none 110 111 !Arguments ------------------------------------ 112 !scalars 113 integer,intent(in) :: icg,icg1,idir,ikpt,ipert,isppol,istwf_k 114 integer,intent(in) :: mkmem,mk1mem,mpert,mpw,mpw1,nsppol 115 integer,intent(inout) :: nband_k,npw1_k,npw_k 116 real(dp),intent(in) :: wtk_k 117 type(MPI_type),intent(in) :: mpi_enreg 118 type(dataset_type),intent(in) :: dtset 119 type(gs_hamiltonian_type),intent(inout) :: gs_hamkq 120 type(pseudopotential_type),intent(in) :: psps 121 !arrays 122 integer,intent(in) :: ddkfil(3),kg1_k(3,npw1_k) 123 integer,intent(in) :: kg_k(3,npw_k) 124 real(dp),intent(in) :: cg(2,mpw*dtset%nspinor*dtset%mband*mkmem*nsppol) 125 real(dp),intent(in) :: cg1(2,mpw1*dtset%nspinor*dtset%mband*mk1mem*nsppol) 126 real(dp),intent(in) :: eig_k(dtset%mband*nsppol),kpt(3),kpq(3),occ_k(nband_k),rmet(3,3) 127 real(dp),intent(in) :: ylm(npw_k,psps%mpsang*psps%mpsang*psps%useylm) 128 real(dp),intent(in) :: ylm1(npw1_k,psps%mpsang*psps%mpsang*psps%useylm) 129 real(dp),intent(inout) :: eig1_k(2*nsppol*dtset%mband**2) 130 real(dp),intent(out) :: d2bbb_k(2,3,dtset%mband,dtset%mband*dtset%prtbbb) 131 real(dp),intent(inout) :: d2nl_k(2,3,mpert) 132 type(wfk_t),intent(inout) :: ddks(3) 133 134 !Local variables------------------------------- 135 !scalars 136 integer :: berryopt,dimffnl,dimffnl1,dimph3d 137 integer :: iband,ider,idir1,ipert1,ipw,jband,nband_kocc,nkpg,nkpg1 !ierr,ii 138 integer :: npw_disk,nsp,optlocal,optnl,opt_gvnl1,sij_opt,tim_getgh1c,usevnl 139 logical :: ddk 140 real(dp) :: aa,dot1i,dot1r,dot2i,dot2r,dot_ndiagi,dot_ndiagr,doti,dotr,lambda 141 character(len=500) :: msg 142 type(rf_hamiltonian_type) :: rf_hamkq 143 !arrays 144 integer :: ik_ddks(3) 145 real(dp) :: dum_grad_berry(1,1),dum_gvnl1(1,1),dum_gs1(1,1),dum_ylmgr(1,3,1),tsec(2) 146 real(dp),allocatable :: cg_k(:,:),cwave0(:,:),cwavef(:,:),cwavef_da(:,:) 147 real(dp),allocatable :: cwavef_db(:,:),dkinpw(:),eig2_k(:),ffnl1(:,:,:,:),ffnlk(:,:,:,:) 148 real(dp),allocatable :: gvnl1(:,:),kinpw1(:),kpg1_k(:,:),kpg_k(:,:),ph3d(:,:,:) 149 type(pawcprj_type),allocatable :: dum_cwaveprj(:,:) 150 151 ! ********************************************************************* 152 153 DBG_ENTER("COLL") 154 155 !Not valid for PAW 156 if (psps%usepaw==1) then 157 msg=' This routine cannot be used for PAW (use pawnst3 instead) !' 158 MSG_BUG(msg) 159 end if 160 161 !Keep track of total time spent in dfpt_nstwf 162 call timab(102,1,tsec) 163 tim_getgh1c=2 164 165 !Miscelaneous inits 166 ABI_DATATYPE_ALLOCATE(dum_cwaveprj,(0,0)) 167 ddk=(ipert==dtset%natom+1.or.ipert==dtset%natom+10.or.ipert==dtset%natom+11) 168 169 !Additional allocations 170 if (.not.ddk) then 171 ABI_ALLOCATE(dkinpw,(npw_k)) 172 ABI_ALLOCATE(kinpw1,(npw1_k)) 173 kinpw1=zero;dkinpw=zero 174 else 175 ABI_ALLOCATE(dkinpw,(0)) 176 ABI_ALLOCATE(kinpw1,(0)) 177 end if 178 ABI_ALLOCATE(gvnl1,(2,npw1_k*dtset%nspinor)) 179 ABI_ALLOCATE(eig2_k,(2*nsppol*dtset%mband**2)) 180 ABI_ALLOCATE(cwave0,(2,npw_k*dtset%nspinor)) 181 ABI_ALLOCATE(cwavef,(2,npw1_k*dtset%nspinor)) 182 183 !Compute (k+G) vectors 184 nkpg=0;if (.not.ddk) nkpg=3*gs_hamkq%nloalg(3) 185 ABI_ALLOCATE(kpg_k,(npw_k,nkpg)) 186 if (nkpg>0) then 187 call mkkpg(kg_k,kpg_k,kpt,nkpg,npw_k) 188 end if 189 190 !Compute (k+q+G) vectors 191 nkpg1=0;if (.not.ddk) nkpg1=3*gs_hamkq%nloalg(3) 192 ABI_ALLOCATE(kpg1_k,(npw1_k,nkpg1)) 193 if (nkpg1>0) then 194 call mkkpg(kg1_k,kpg1_k,kpq,nkpg1,npw1_k) 195 end if 196 197 !Compute nonlocal form factors ffnl at (k+G) 198 dimffnl=0;if (.not.ddk) dimffnl=1 199 ABI_ALLOCATE(ffnlk,(npw_k,dimffnl,psps%lmnmax,psps%ntypat)) 200 if (.not.ddk) then 201 ider=0 202 call mkffnl(psps%dimekb,dimffnl,psps%ekb,ffnlk,psps%ffspl,gs_hamkq%gmet,& 203 & gs_hamkq%gprimd,ider,ider,psps%indlmn,kg_k,kpg_k,kpt,psps%lmnmax,& 204 & psps%lnmax,psps%mpsang,psps%mqgrid_ff,nkpg,npw_k,psps%ntypat,psps%pspso,& 205 & psps%qgrid_ff,rmet,psps%usepaw,psps%useylm,ylm,dum_ylmgr) 206 end if 207 208 !Compute nonlocal form factors ffnl1 at (k+q+G) 209 dimffnl1=0;if (.not.ddk) dimffnl1=1 210 ABI_ALLOCATE(ffnl1,(npw1_k,dimffnl1,psps%lmnmax,psps%ntypat)) 211 if (.not.ddk) then 212 ider=0 213 call mkffnl(psps%dimekb,dimffnl1,psps%ekb,ffnl1,psps%ffspl,gs_hamkq%gmet,& 214 & gs_hamkq%gprimd,ider,ider,psps%indlmn,kg1_k,kpg1_k,kpq,& 215 & psps%lmnmax,psps%lnmax,psps%mpsang,psps%mqgrid_ff,nkpg1,npw1_k,psps%ntypat,& 216 & psps%pspso,psps%qgrid_ff,rmet,psps%usepaw,psps%useylm,ylm1,dum_ylmgr) 217 end if 218 219 !Load k-dependent part in the Hamiltonian datastructure 220 call load_k_hamiltonian(gs_hamkq,kpt_k=kpt,npw_k=npw_k,istwf_k=istwf_k,& 221 & kg_k=kg_k,kpg_k=kpg_k,ffnl_k=ffnlk) 222 223 !Load k+q-dependent part in the Hamiltonian datastructure 224 ! Note: istwf_k is imposed to 1 for RF calculations (should use istwf_kq instead) 225 dimph3d=0;if (.not.ddk) dimph3d=gs_hamkq%matblk 226 ABI_ALLOCATE(ph3d,(2,npw1_k,dimph3d)) 227 call load_kprime_hamiltonian(gs_hamkq,kpt_kp=kpq,npw_kp=npw1_k,istwf_kp=istwf_k,& 228 & kinpw_kp=kinpw1,kg_kp=kg1_k,kpg_kp=kpg1_k,ffnl_kp=ffnl1,& 229 & ph3d_kp=ph3d,compute_ph3d=(.not.ddk)) 230 231 !Load k-dependent part in the 1st-order Hamiltonian datastructure 232 call load_k_rf_hamiltonian(rf_hamkq,npw_k=npw_k,dkinpw_k=dkinpw) 233 234 !Take care of the npw and kg records 235 !NOTE : one should be able to modify the rwwf routine to take care 236 !of the band parallelism, which is not the case yet ... 237 ik_ddks = 0 238 do idir1=1,3 239 if (ddkfil(idir1)/=0)then 240 ! Read npw record 241 nsp=dtset%nspinor 242 ik_ddks(idir1) = wfk_findk(ddks(idir1), kpt) 243 ABI_CHECK(ik_ddks(idir1) /= -1, "Cannot find kpt") 244 npw_disk = ddks(idir1)%hdr%npwarr(ik_ddks(idir1)) 245 if (npw_k /= npw_disk) then 246 write(unit=msg,fmt='(a,i3,a,i5,a,i3,a,a,i5,a,a,i5)')& 247 & 'For isppol = ',isppol,', ikpt = ',ikpt,' and idir = ',idir,ch10,& 248 & 'the number of plane waves in the ddk file is equal to', npw_disk,ch10,& 249 & 'while it should be ',npw_k 250 MSG_BUG(msg) 251 end if 252 end if 253 end do 254 255 if (ipert==dtset%natom+1) then 256 nband_kocc = 0 257 do iband = 1,nband_k 258 if (abs(occ_k(iband)) > tol8) nband_kocc = nband_kocc + 1 259 nband_kocc = max (nband_kocc, 1) 260 end do 261 end if 262 263 if(dtset%prtbbb==1)then 264 ABI_ALLOCATE(cwavef_da,(2,npw1_k*dtset%nspinor)) 265 ABI_ALLOCATE(cwavef_db,(2,npw1_k*dtset%nspinor)) 266 ABI_ALLOCATE(cg_k,(2,npw_k*dtset%nspinor*nband_k)) 267 if ((ipert == dtset%natom + 1).or.(ipert <= dtset%natom).or. & 268 & (ipert == dtset%natom + 2).or.(ipert == dtset%natom + 5)) then 269 cg_k(:,:) = cg(:,1+icg:icg+nband_k*npw_k*dtset%nspinor) 270 end if 271 d2bbb_k(:,:,:,:) = zero 272 end if 273 274 !Loop over bands 275 do iband=1,nband_k 276 277 if(mpi_enreg%proc_distrb(ikpt,iband,isppol) /= mpi_enreg%me_kpt) cycle 278 279 ! Read ground-state wavefunctions 280 if (dtset%prtbbb==0 .or. ipert==dtset%natom+2) then 281 cwave0(:,:)=cg(:,1+(iband-1)*npw_k*dtset%nspinor+icg:iband*npw_k*dtset%nspinor+icg) 282 else ! prtbbb==1 and ipert<=natom , already in cg_k 283 cwave0(:,:)=cg_k(:,1+(iband-1)*npw_k*dtset%nspinor:iband*npw_k*dtset%nspinor) 284 end if 285 286 ! Get first-order wavefunctions 287 cwavef(:,:)=cg1(:,1+(iband-1)*npw1_k*dtset%nspinor+icg1:iband*npw1_k*dtset%nspinor+icg1) 288 289 ! In case non ddk perturbation 290 if (ipert /= dtset%natom + 1) then 291 292 do ipert1=1,mpert 293 294 if( ipert1<=dtset%natom .or. ipert1==dtset%natom+2 )then 295 296 ! Initialize data for NL 1st-order hamiltonian 297 call init_rf_hamiltonian(1,gs_hamkq,ipert1,rf_hamkq) 298 299 if (((ipert <= dtset%natom).or.(ipert == dtset%natom + 2)) & 300 & .and.(ipert1 == dtset%natom+2).and. dtset%prtbbb==1) then 301 call gaugetransfo(cg_k,cwavef,cwavef_db,eig_k,eig1_k,iband,nband_k, & 302 & dtset%mband,npw_k,npw1_k,dtset%nspinor,nsppol,occ_k) 303 cwavef(:,:) = cwavef_db(:,:) 304 end if 305 306 ! Define the direction along which to move the atom : 307 ! the polarisation (ipert1,idir1) is refered as j1. 308 do idir1=1,3 309 if (ipert1<=dtset%natom.or.(ipert1==dtset%natom+2.and.ddkfil(idir1)/=0)) then 310 311 ! Get |Vnon-locj^(1)|u0> : 312 ! First-order non-local, applied to zero-order wavefunction 313 ! This routine gives MINUS the non-local contribution 314 315 ! ==== Atomic displ. perturbation 316 if( ipert1<=dtset%natom )then 317 lambda=eig_k((isppol-1)*nband_k+iband) 318 berryopt=1;optlocal=0;optnl=1;usevnl=0;opt_gvnl1=0;sij_opt=0 319 call getgh1c(berryopt,cwave0,dum_cwaveprj,gvnl1,dum_grad_berry,& 320 & dum_gs1,gs_hamkq,dum_gvnl1,idir1,ipert1,lambda,mpi_enreg,optlocal,& 321 & optnl,opt_gvnl1,rf_hamkq,sij_opt,tim_getgh1c,usevnl) 322 323 ! ==== Electric field perturbation 324 else if( ipert1==dtset%natom+2 )then 325 ! TODO: Several tests fail here ifdef HAVE_MPI_IO_DEFAULT 326 ! The problem is somehow related to the use of MPI-IO file views!. 327 call wfk_read_bks(ddks(idir1), iband, ik_ddks(idir1), isppol, xmpio_single, cg_bks=gvnl1, & 328 eig1_bks=eig2_k(1+(iband-1)*2*nband_k:)) 329 !eig1_bks=eig2_k(1+(iband-1)*2*nband_k:2*iband*nband_k)) 330 !write(777,*)"eig2_k, gvnl1 for band: ",iband,", ikpt: ",ikpt 331 !do ii=1,2*nband_k 332 ! write(777,*)eig2_k(ii+(iband-1)) 333 !end do 334 !write(777,*)gvnl1 335 336 ! In case of band-by-band, 337 ! construct the first-order wavefunctions in the diagonal gauge 338 if (((ipert <= dtset%natom).or.(ipert == dtset%natom + 2)).and.(dtset%prtbbb==1)) then 339 call gaugetransfo(cg_k,gvnl1,cwavef_da,eig_k,eig2_k,iband,nband_k, & 340 & dtset%mband,npw_k,npw1_k,dtset%nspinor,nsppol,occ_k) 341 gvnl1(:,:) = cwavef_da(:,:) 342 end if 343 ! Multiplication by -i 344 do ipw=1,npw1_k*dtset%nspinor 345 aa=gvnl1(1,ipw) 346 gvnl1(1,ipw)=gvnl1(2,ipw) 347 gvnl1(2,ipw)=-aa 348 end do 349 end if 350 351 ! MVeithen 021212 : 352 ! 1) Case ipert1 = natom + 2 and ipert = natom + 2: 353 ! the second derivative of the energy with respect to an electric 354 ! field is computed from Eq. (38) of X. Gonze, PRB 55 ,10355 (1997). 355 ! The evaluation of this formula needs the operator $i \frac{d}{dk}. 356 ! 2) Case ipert1 = natom + 2 and ipert < natom: 357 ! the computation of the Born effective charge tensor uses 358 ! the operator $-i \frac{d}{dk}. 359 if (ipert==dtset%natom+2) gvnl1(:,:) = -gvnl1(:,:) 360 361 ! <G|Vnl1|Cnk> is contained in gvnl1 362 ! construct the matrix element (<uj2|vj1|u0>)complex conjug and add it to the 2nd-order matrix 363 call dotprod_g(dotr,doti,istwf_k,npw1_k*dtset%nspinor,2,cwavef,gvnl1,mpi_enreg%me_g0,mpi_enreg%comm_spinorfft) 364 d2nl_k(1,idir1,ipert1)=d2nl_k(1,idir1,ipert1)+wtk_k*occ_k(iband)*two*dotr 365 d2nl_k(2,idir1,ipert1)=d2nl_k(2,idir1,ipert1)-wtk_k*occ_k(iband)*two*doti 366 367 ! Band by band decomposition of the Born effective charges 368 ! calculated from a phonon perturbation 369 if(dtset%prtbbb==1)then 370 d2bbb_k(1,idir1,iband,iband) = wtk_k*occ_k(iband)*two*dotr 371 d2bbb_k(2,idir1,iband,iband) = -one*wtk_k*occ_k(iband)*two*doti 372 end if 373 374 end if 375 end do 376 377 call destroy_rf_hamiltonian(rf_hamkq) 378 end if 379 end do 380 end if ! ipert /= natom +1 381 382 ! Compute the localization tensor 383 384 if (ipert==dtset%natom+1) then 385 386 ipert1=dtset%natom+1 387 if(dtset%prtbbb==1)then 388 call gaugetransfo(cg_k,cwavef,cwavef_db,eig_k,eig1_k,iband,nband_k, & 389 & dtset%mband,npw_k,npw1_k,dtset%nspinor,nsppol,occ_k) 390 cwavef(:,:) = cwavef_db(:,:) 391 end if 392 393 do idir1 = 1,3 394 eig2_k(:) = zero 395 gvnl1(:,:) = zero 396 if (idir == idir1) then 397 gvnl1(:,:) = cwavef(:,:) 398 eig2_k(:) = eig1_k(:) 399 else 400 if (ddkfil(idir1) /= 0) then 401 call wfk_read_bks(ddks(idir1), iband, ik_ddks(idir1), isppol, xmpio_single, cg_bks=gvnl1, & 402 eig1_bks=eig2_k(1+(iband-1)*2*nband_k:)) 403 !eig1_bks=eig2_k(1+(iband-1)*2*nband_k:2*iband*nband_k)) 404 !write(778,*)"eig2_k, gvnl1 for band: ",iband,", ikpt: ",ikpt 405 !do ii=1,2*nband_k 406 ! write(778,*)eig2_k(ii+(iband-1)) 407 !end do 408 !write(778,*)gvnl1 409 410 if(dtset%prtbbb==1)then 411 call gaugetransfo(cg_k,gvnl1,cwavef_da,eig_k,eig2_k,iband,nband_k, & 412 & dtset%mband,npw_k,npw1_k,dtset%nspinor,nsppol,occ_k) 413 gvnl1(:,:) = cwavef_da(:,:) 414 end if 415 416 end if !ddkfil(idir1) 417 end if !idir == idir1 418 419 ! <G|du/dqa> is contained in gvnl1 and <G|du/dqb> in cwavef 420 ! construct the matrix elements <du/dqa|du/dqb> -> dot 421 ! <u|du/dqa> -> dot1 422 ! <du/dqb|u> -> dot2 423 ! and add them to the 2nd-order matrix 424 425 call dotprod_g(dotr,doti,istwf_k,npw1_k*dtset%nspinor,2,gvnl1,cwavef,mpi_enreg%me_g0,mpi_enreg%comm_spinorfft) 426 d2nl_k(1,idir1,ipert1)=d2nl_k(1,idir1,ipert1)+wtk_k*occ_k(iband)*dotr/(nband_kocc*two) 427 d2nl_k(2,idir1,ipert1)=d2nl_k(2,idir1,ipert1)+wtk_k*occ_k(iband)*doti/(nband_kocc*two) 428 429 430 ! XG 020216 : Marek, could you check the next forty lines 431 ! In the parallel gauge, dot1 and dot2 vanishes 432 if(dtset%prtbbb==1)then 433 d2bbb_k(1,idir1,iband,iband)=d2bbb_k(1,idir1,iband,iband)+dotr 434 d2bbb_k(2,idir1,iband,iband)=d2bbb_k(2,idir1,iband,iband)+doti 435 dot_ndiagr=zero ; dot_ndiagi=zero 436 do jband = 1,nband_k !compute dot1 and dot2 437 if (abs(occ_k(jband)) > tol8) then 438 dot1r=zero ; dot1i=zero 439 dot2r=zero ; dot2i=zero 440 cwave0(:,:)=cg_k(:,1+(jband-1)*npw_k*dtset%nspinor:jband*npw_k*dtset%nspinor) 441 442 call dotprod_g(dot1r,dot1i,istwf_k,npw1_k*dtset%nspinor,2,cwave0,gvnl1,mpi_enreg%me_g0,mpi_enreg%comm_spinorfft) 443 call dotprod_g(dot2r,dot2i,istwf_k,npw1_k*dtset%nspinor,2,cwavef,cwave0,mpi_enreg%me_g0,mpi_enreg%comm_spinorfft) 444 445 dot_ndiagr = dot_ndiagr + dot1r*dot2r - dot1i*dot2i 446 dot_ndiagi = dot_ndiagi + dot1r*dot2i + dot1i*dot2r 447 d2bbb_k(1,idir1,iband,jband) = d2bbb_k(1,idir1,iband,jband) - & 448 & (dot1r*dot2r - dot1i*dot2i) 449 d2bbb_k(2,idir1,iband,jband) = d2bbb_k(2,idir1,iband,jband) - & 450 & (dot1r*dot2i + dot1i*dot2r) 451 end if ! occ_k 452 end do !jband 453 d2bbb_k(:,idir1,iband,:)=d2bbb_k(:,idir1,iband,:)*wtk_k*occ_k(iband)*half 454 d2nl_k(1,idir1,ipert1)= & 455 & d2nl_k(1,idir1,ipert1)-wtk_k*occ_k(iband)*dot_ndiagr/(nband_kocc*two) 456 d2nl_k(2,idir1,ipert1)=& 457 & d2nl_k(2,idir1,ipert1)-wtk_k*occ_k(iband)*dot_ndiagi/(nband_kocc*two) 458 end if ! prtbbb==1 459 460 end do ! idir1 461 end if ! Compute localization tensor, ipert=natom+1 462 463 ! End loop over bands 464 end do 465 466 !Final deallocations 467 ABI_DEALLOCATE(cwave0) 468 ABI_DEALLOCATE(cwavef) 469 ABI_DEALLOCATE(eig2_k) 470 ABI_DEALLOCATE(gvnl1) 471 ABI_DEALLOCATE(ffnlk) 472 ABI_DEALLOCATE(ffnl1) 473 ABI_DEALLOCATE(dkinpw) 474 ABI_DEALLOCATE(kinpw1) 475 ABI_DEALLOCATE(kpg_k) 476 ABI_DEALLOCATE(kpg1_k) 477 ABI_DEALLOCATE(ph3d) 478 ABI_DATATYPE_DEALLOCATE(dum_cwaveprj) 479 if(dtset%prtbbb==1) then 480 ABI_DEALLOCATE(cg_k) 481 ABI_DEALLOCATE(cwavef_da) 482 ABI_DEALLOCATE(cwavef_db) 483 end if 484 485 call timab(102,2,tsec) 486 487 DBG_EXIT("COLL") 488 489 end subroutine dfpt_nstwf