TABLE OF CONTENTS
ABINIT/eig2stern [ Functions ]
NAME
eig2stern
FUNCTION
This routine calculates the second-order eigenvalues. The output eig2nkq is this quantity for the input k points.
COPYRIGHT
Copyright (C) 1999-2018 ABINIT group (PB, XG) 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 .
INPUTS
bdeigrf = number of bands for which to calculate the second-order eigenvalues.
clflg(3,mpert)= array on calculated perturbations for eig2rf.
dim_eig2nkq = 1 if eig2nkq is to be computed.
cg1_pert(2,mpw1*nspinor*mband*mk1mem*nsppol,3,mpert) = first-order wf in G
space for each perturbation. The wavefunction is orthogonal to the
active space.
gh0c1_pert(2,mpw1*nspinor*mband*mk1mem*nsppol,3,mpert) = matrix containing the
vector: <G|H(0)|psi(1)>, for each perturbation.
gh1c_pert(2,mpw1*nspinor*mband*mk1mem*nsppol,3,mpert)) = matrix containing the
vector: <G|H(1)|n,k>, for each perturbation. The wavefunction is
orthogonal to the active space.
eigbrd(2,mband*nsppol,nkpt,3,natom,3,natom) = broadening factors for the
electronic eigenvalues (optional).
eigen0(nkpt_rbz*mband*nsppol) = 0-order eigenvalues at all K-points:
<k,n'|H(0)|k,n'> (hartree).
eigenq(nkpt_rbz*mband*nsppol) = 0-order eigenvalues at all shifted K-points:
<k+Q,n'|H(0)|k+Q,n'> (hartree).
eigen1(nkpt_rbz*2*nsppol*mband**2,3,mpert) = matrix of first-order:
<k+Q,n'|H(1)|k,n> (hartree) (calculated in dfpt_cgwf).
eig2nkq(2,mband*nsppol,nkpt,3,natom,3,natom*dim_eig2nkq) = second derivatives of
the electronic eigenvalues.
elph2_imagden = imaginary part of the denominator of the sum-over-state expression
for the electronic eigenenergy shift due to second-order electron-phonon
interation.
ieig2rf = integer for calculation type.
indsym(4,nsym,natom) = indirect indexing array for atom labels
(not used yet, but will be used with symmetries).
istwfk_pert(nkpt_rbz,3,mpert) = integer for choice of storage of wavefunction at
each k point for each perturbation.
mband = maximum number of bands.
mk1mem = maximum number of k points which can fit in memory (RF data);
0 if use disk.
mpert = maximum number of perturbations.
natom = number of atoms in the unit cell.
npert = number of phonon perturbations, without taking into account directions:
natom.
nsym = number of symmetries (not used yet).
mpi_enreg = informations about MPI parallelization.
mpw1 = maximum number of planewaves used to represent first-order wavefunctions.
nkpt_rbz = number of k-points for each perturbation.
npwar1(nkpt_rbz,mpert) = number of planewaves at k-point for first-order.
nspinor = number of spinorial components of the wavefunctions.
nsppol = 1 for unpolarized, 2 for spin-polarized.
occ(mband*nkpt*nsppol)=occup number for each band (often 2) at each k point
smdelta = integer controling the calculation of electron lifetimes.
symq(4,2,nsym) = 1 if symmetry preserves present qpoint. From littlegroup_q (not used yet).
symrec(3,3,nsym) = 3x3 matrices of the group symmetries (reciprocal space)
(not used yet).
symrel(3,3,nsym) = array containing the symmetries in real space (not used yet).
timrev = 1 if time-reversal preserves the q wavevector; 0 otherwise
(not in use yet).
dtset = OPTIONAL, dataset structure containing the input variable of the
calculation. This is required to use the k-interpolation routine.
eigenq_fine(mband_fine,mkpt_fine,nsppol_fine) = OPTIONAL, 0-order eigenvalues
at all shifted K-points: <k+Q,n'|H(0)|k+Q,n'> (hartree) of the
fine grid. This information is read from the WF dense k-grid file.
hdr_fine = OPTIONAL, header of the WF file of the fine k-point grid. This
variable is required for the k-interpolation routine.
hdr0 = OPTIONAL, header of the GS WF file of the corse k-point grid. This
variable is required for the k-interpolation routine.
OUTPUT
eig2nkq(2,mband*nsppol,nkpt_rbz,3,npert,3,npert)= diagonal part of the
second-order eigenvalues: E^{(2),diag}_{k,q,j}.
eigbrd(2,mband*nsppol,nkpt_rbz,3,npert,3,npert)= OPTIONAL, array containing the
electron lifetimes.
PARENTS
dfpt_looppert
CHILDREN
distrb2,dotprod_g,kptfine_av,smeared_delta,timab,wrtout,xmpi_sum
SOURCE
93 #if defined HAVE_CONFIG_H 94 #include "config.h" 95 #endif 96 97 #include "abi_common.h" 98 99 subroutine eig2stern(occ,bdeigrf,clflg,cg1_pert,dim_eig2nkq,dim_eig2rf,eigen0,eigenq,& 100 & eigen1,eig2nkq,elph2_imagden,esmear,gh0c1_pert,gh1c_pert,ieig2rf,istwfk_pert,& 101 & mband,mk1mem,mpert,npert,mpi_enreg,mpw1,nkpt_rbz,npwar1,nspinor,nsppol,smdelta,& 102 & dtset,eigbrd,eigenq_fine,hdr_fine,hdr0) 103 104 use defs_basis 105 use defs_abitypes 106 use m_xmpi 107 use m_errors 108 use m_cgtools 109 use m_profiling_abi 110 111 !This section has been created automatically by the script Abilint (TD). 112 !Do not modify the following lines by hand. 113 #undef ABI_FUNC 114 #define ABI_FUNC 'eig2stern' 115 use interfaces_14_hidewrite 116 use interfaces_18_timing 117 use interfaces_32_util 118 use interfaces_51_manage_mpi 119 use interfaces_72_response, except_this_one => eig2stern 120 !End of the abilint section 121 122 implicit none 123 124 !Arguments ------------------------------------ 125 !scalars 126 integer,intent(in) :: bdeigrf,dim_eig2nkq,dim_eig2rf,ieig2rf,mband,mk1mem,mpert,mpw1,nkpt_rbz 127 integer,intent(in) :: npert,nspinor,nsppol,smdelta 128 real(dp),intent(in) :: elph2_imagden,esmear 129 type(MPI_type),intent(inout) :: mpi_enreg 130 !arrays 131 integer,intent(in) :: clflg(3,mpert) 132 integer,intent(in) :: istwfk_pert(nkpt_rbz,3,mpert) 133 integer,intent(in) :: npwar1(nkpt_rbz,mpert) 134 real(dp),intent(in) :: cg1_pert(2,mpw1*nspinor*mband*mk1mem*nsppol*dim_eig2rf,3,mpert) 135 real(dp),intent(in) :: gh0c1_pert(2,mpw1*nspinor*mband*mk1mem*nsppol*dim_eig2rf,3,mpert) 136 real(dp),intent(in) :: gh1c_pert(2,mpw1*nspinor*mband*mk1mem*nsppol*dim_eig2rf,3,mpert) 137 real(dp),intent(inout) :: eigen0(nkpt_rbz*mband*nsppol) 138 real(dp),intent(in) :: eigen1(nkpt_rbz*2*nsppol*mband**2,3,mpert) 139 real(dp),intent(inout) :: eigenq(nkpt_rbz*mband*nsppol) 140 real(dp),intent(out) :: eig2nkq(2,mband*nsppol,nkpt_rbz,3,npert,3,npert*dim_eig2nkq) 141 real(dp),intent(out),optional :: eigbrd(2,mband*nsppol,nkpt_rbz,3,npert,3,npert) 142 real(dp),intent(in),pointer,optional :: eigenq_fine(:,:,:) 143 real(dp), intent(in) :: occ(mband*nkpt_rbz*nsppol) 144 type(dataset_type), intent(in) :: dtset 145 type(hdr_type),intent(in),optional :: hdr_fine,hdr0 146 147 !Local variables------------------------------- 148 !tolerance for non degenerated levels 149 !scalars 150 integer :: band2tot_index,band_index,bandtot_index,iband,icg2,idir1,idir2 151 integer :: ikpt,ipert1,ipert2,isppol,istwf_k,jband,npw1_k,nkpt_sub,ikpt2 152 !integer :: ipw 153 integer :: master,me,spaceworld,ierr 154 !real(dp),parameter :: etol=1.0d-3 155 real(dp),parameter :: etol=1.0d-6 156 !real(dp),parameter :: etol=zero 157 real(dp) :: ar,ai,deltae,den,dot2i,dot2r,dot3i,dot3r,doti,dotr,eig1_i1,eig1_i2 158 real(dp) :: eig1_r1,eig1_r2,eig2_diai,den_av 159 real(dp) :: wgt_int 160 real(dp) :: eig2_diar,eigbrd_i,eigbrd_r 161 character(len=500) :: message 162 character(len=500) :: msg 163 !DBSP 164 ! character(len=300000) :: message2 165 !END 166 logical :: test_do_band 167 !arrays 168 integer, allocatable :: nband_rbz(:),icg2_rbz(:,:) 169 integer,pointer :: kpt_fine_sub(:) 170 real(dp) :: tsec(2) 171 real(dp),allocatable :: cwavef(:,:),cwavef2(:,:),center(:),eigen0tmp(:),eigenqtmp(:) 172 real(dp) :: eigen(mband*nsppol),eigen_prime(mband*nsppol) 173 real(dp),allocatable :: gh(:,:),gh1(:,:),ghc(:,:) 174 real(dp),allocatable :: smdfun(:,:) 175 real(dp),pointer :: wgt_sub(:) 176 177 ! ********************************************************************* 178 179 !Init parallelism 180 master =0 181 spaceworld=mpi_enreg%comm_cell 182 me=mpi_enreg%me_kpt 183 !DEBUG 184 !write(std_out,*)' eig2stern : enter ' 185 !write(std_out,*)' mpw1=',mpw1 186 !write(std_out,*)' mband=',mband 187 !write(std_out,*)' nsppol=',nsppol 188 !write(std_out,*)' nkpt_rbz=',nkpt_rbz 189 !write(std_out,*)' npert=',npert 190 !ENDDEBUG 191 192 !Init interpolation method 193 if(present(eigenq_fine))then 194 ABI_ALLOCATE(center,(3)) 195 end if 196 197 call timab(148,1,tsec) 198 199 if(nsppol==2)then 200 message = 'nsppol=2 is still under development. Be careful when using it ...' 201 MSG_COMMENT(message) 202 end if 203 204 band2tot_index =0 205 bandtot_index=0 206 band_index=0 207 208 !Add scissor shift to eigenenergies 209 if (dtset%dfpt_sciss > tol6 ) then 210 write(msg,'(a,f7.3,2a)')& 211 & ' A scissor operator of ',dtset%dfpt_sciss*Ha_eV,' [eV] has been applied to the eigenenergies',ch10 212 call wrtout(std_out,msg,'COLL') 213 call wrtout(ab_out,msg,'COLL') 214 ABI_ALLOCATE(eigen0tmp,(nkpt_rbz*mband*nsppol)) 215 ABI_ALLOCATE(eigenqtmp,(nkpt_rbz*mband*nsppol)) 216 eigen0tmp = eigen0(:) 217 eigenqtmp = eigenq(:) 218 eigen0 = zero 219 eigenq = zero 220 end if 221 222 if(ieig2rf > 0) then 223 eig2nkq(:,:,:,:,:,:,:) = zero 224 end if 225 if(present(eigbrd))then 226 eigbrd(:,:,:,:,:,:,:) = zero 227 end if 228 229 if(xmpi_paral==1) then 230 ABI_ALLOCATE(mpi_enreg%proc_distrb,(nkpt_rbz,mband,nsppol)) 231 ABI_ALLOCATE(nband_rbz,(nkpt_rbz*nsppol)) 232 if (allocated(mpi_enreg%my_kpttab)) then 233 ABI_DEALLOCATE(mpi_enreg%my_kpttab) 234 end if 235 ABI_ALLOCATE(mpi_enreg%my_kpttab,(nkpt_rbz)) 236 ! Assume the number of bands is the same for all k points. 237 nband_rbz(:)=mband 238 call distrb2(mband,nband_rbz,nkpt_rbz,mpi_enreg%nproc_cell,nsppol,mpi_enreg) 239 end if 240 241 icg2=0 242 ipert1=1 ! Suppose that the situation is the same for all perturbations 243 ABI_ALLOCATE(icg2_rbz,(nkpt_rbz,nsppol)) 244 do isppol=1,nsppol 245 do ikpt=1,nkpt_rbz 246 icg2_rbz(ikpt,isppol)=icg2 247 if(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,mband,isppol,me)) cycle 248 icg2 = icg2 + npwar1(ikpt,ipert1)*nspinor*mband 249 end do 250 end do 251 252 do isppol=1,nsppol 253 do ikpt =1,nkpt_rbz 254 255 if(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,mband,isppol,me)) then 256 band2tot_index = band2tot_index + 2*mband**2 257 bandtot_index = bandtot_index + mband 258 cycle 259 end if 260 261 if(present(eigenq_fine))then 262 write(std_out,*) 'Start of the energy denominator interpolation method.' 263 nkpt_sub = 0 264 ! center is the k+q point around which we will average the kpt_fine 265 center = hdr0%kptns(:,ikpt)+ dtset%qptn(:) 266 267 call kptfine_av(center,dtset%qptrlatt,hdr_fine%kptns,hdr_fine%nkpt,& 268 & kpt_fine_sub,nkpt_sub,wgt_sub) 269 write(std_out,'(a,3f8.4,a,i3)') 'Number of k-points of the fine grid & 270 & around the k+Q point ',center,' is:',nkpt_sub 271 write(std_out,'(a,f10.5)') 'The sum of the weights of the k-points is: ',SUM(wgt_sub) 272 end if 273 274 ! Add scissor shift to eigenenergies 275 if (dtset%dfpt_sciss > tol6 ) then 276 do iband=1,mband 277 if (occ(iband+bandtot_index) < tol6) then 278 eigen0(iband+bandtot_index) = eigen0tmp(iband+bandtot_index) + dtset%dfpt_sciss 279 eigenq(iband+bandtot_index) = eigenqtmp(iband+bandtot_index) + dtset%dfpt_sciss 280 else 281 eigen0(iband+bandtot_index) = eigen0tmp(iband+bandtot_index) 282 eigenq(iband+bandtot_index) = eigenqtmp(iband+bandtot_index) 283 end if 284 end do 285 end if 286 287 288 if(smdelta >0) then !broadening 289 if(.not.allocated(smdfun)) then 290 ABI_ALLOCATE(smdfun,(mband,mband)) 291 end if 292 smdfun(:,:) = zero 293 do iband=1,mband 294 eigen(iband) = eigen0(iband+bandtot_index) 295 eigen_prime(iband) =eigenq(iband+bandtot_index) 296 end do 297 if(esmear>tol6) then 298 call smeared_delta(eigen,eigen_prime,esmear,mband,smdelta,smdfun) 299 end if 300 end if 301 icg2=icg2_rbz(ikpt,isppol) 302 303 ipert1=1 ! Suppose all perturbations lead to the same number of planewaves 304 npw1_k = npwar1(ikpt,ipert1) 305 ABI_ALLOCATE(cwavef,(2,npw1_k*nspinor)) 306 ABI_ALLOCATE(cwavef2,(2,npw1_k*nspinor)) 307 ABI_ALLOCATE(gh,(2,npw1_k*nspinor)) 308 ABI_ALLOCATE(gh1,(2,npw1_k*nspinor)) 309 ABI_ALLOCATE(ghc,(2,npw1_k*nspinor)) 310 311 do iband=1,bdeigrf 312 313 ! If the k point and band belong to me, compute the contribution 314 test_do_band=.true. 315 if(mpi_enreg%proc_distrb(ikpt,iband,isppol)/=me)test_do_band=.false. 316 317 if(test_do_band)then 318 319 do ipert1=1,npert 320 321 do idir1=1,3 322 if(clflg(idir1,ipert1)==0)cycle 323 istwf_k = istwfk_pert(ikpt,idir1,ipert1) 324 325 do ipert2=1,npert 326 do idir2=1,3 327 if(clflg(idir2,ipert2)==0)cycle 328 329 eig2_diar = zero ; eig2_diai = zero ; eigbrd_r = zero ; eigbrd_i = zero 330 331 do jband=1,mband 332 eig1_r1 = eigen1(2*jband-1+(iband-1)*2*mband+band2tot_index,idir1,ipert1) 333 eig1_r2 = eigen1(2*jband-1+(iband-1)*2*mband+band2tot_index,idir2,ipert2) 334 eig1_i1 = eigen1(2*jband+(iband-1)*2*mband+band2tot_index,idir1,ipert1) 335 eig1_i2 = - eigen1(2*jband+(iband-1)*2*mband+band2tot_index,idir2,ipert2) !the negative sign is from the CC 336 ! If no interpolation, fallback on to the previous 337 ! implementation 338 if(.not. present(eigenq_fine))then 339 deltae=eigenq(jband+bandtot_index)-eigen0(iband+bandtot_index) 340 end if 341 ar=eig1_r1*eig1_r2-eig1_i1*eig1_i2 342 ai=eig1_r1*eig1_i2+eig1_i1*eig1_r2 343 344 ! Sum over all active space to retrieve the diagonal gauge 345 if(ieig2rf == 1 .or. ieig2rf ==2 ) then 346 ! if(abs(deltae)>etol) then ! This is commented because 347 ! there is no problem with divergencies with elph2_imag != 0 348 if( present(eigenq_fine))then 349 den_av = zero 350 wgt_int = zero 351 do ikpt2=1,nkpt_sub 352 deltae=eigenq_fine(jband,kpt_fine_sub(ikpt2),1)& 353 & -eigen0(iband+bandtot_index) 354 den_av = den_av-(wgt_sub(ikpt2)*deltae)/(deltae**2+elph2_imagden**2) 355 wgt_int = wgt_int+wgt_sub(ikpt2) 356 end do 357 den = den_av/wgt_int 358 else 359 if(abs(elph2_imagden) < etol) then 360 if(abs(deltae)>etol) then 361 den=-one/(deltae**2+elph2_imagden**2) 362 else 363 den= zero 364 end if 365 else 366 den=-one/(deltae**2+elph2_imagden**2) 367 end if 368 end if 369 370 ! The following should be the most general implementation of the presence of elph2_imagden 371 ! eig2_diar=eig2_diar+(ar*deltae+ai*elph2_imagden)*den 372 ! eig2_diai=eig2_diai+(ai*deltae-ar*elph2_imagden)*den 373 ! This gives back the implementation without elph2_imagden 374 ! eig2_diar=eig2_diar+ar*deltae*den 375 ! eig2_diai=eig2_diai+ai*deltae*den 376 ! This is what Samuel had implemented 377 ! eig2_diar=eig2_diar+ar*deltae*den 378 ! eig2_diai=eig2_diai+ai*elph2_imagden*den 379 ! Other possibility : throw away the broadening part, that is actually treated separately. 380 if( present(eigenq_fine))then 381 eig2_diar=eig2_diar+ar*den 382 eig2_diai=eig2_diai+ai*den 383 else 384 eig2_diar=eig2_diar+ar*deltae*den 385 eig2_diai=eig2_diai+ai*deltae*den 386 !DBSP 387 ! if (iband+band_index==2 .and. ikpt==1 .and. idir1==1 .and. ipert1==1 .and. idir2==1 .and. ipert2==1) then 388 ! write(message2,*) 'eig2_diar1=',eig2_diar,' ar=',ar,' deltae=',deltae,' den=',den 389 ! call wrtout(std_out,message2,'PERS') 390 ! endif 391 !END 392 393 end if 394 end if ! ieig2rf==1 or 2 395 396 if(present(eigbrd))then 397 if(smdelta >0) then !broadening 398 eigbrd_r = eigbrd_r + ar*smdfun(iband,jband) 399 eigbrd_i = eigbrd_i + ai*smdfun(iband,jband) 400 end if 401 end if 402 403 end do !jband 404 405 ! Add the contribution of non-active bands, if DFPT calculation (= Sternheimer) 406 if(ieig2rf == 1 .or. ieig2rf ==3 .or. ieig2rf ==4 .or. ieig2rf==5 ) then 407 ! if(ieig2rf == 1 ) then 408 409 dotr=zero ; doti=zero 410 dot2r=zero ; dot2i=zero 411 dot3r=zero ; dot3i=zero 412 413 414 cwavef(:,:) = cg1_pert(:,1+(iband-1)*npw1_k*nspinor+icg2:iband*npw1_k*nspinor+icg2,idir2,ipert2) 415 cwavef2(:,:)= cg1_pert(:,1+(iband-1)*npw1_k*nspinor+icg2:iband*npw1_k*nspinor+icg2,idir1,ipert1) 416 gh1(:,:) = gh1c_pert(:,1+(iband-1)*npw1_k*nspinor+icg2:iband*npw1_k*nspinor+icg2,idir1,ipert1) 417 gh(:,:) = gh1c_pert(:,1+(iband-1)*npw1_k*nspinor+icg2:iband*npw1_k*nspinor+icg2,idir2,ipert2) 418 ghc(:,:) = gh0c1_pert(:,1+(iband-1)*npw1_k*nspinor+icg2:iband*npw1_k*nspinor+icg2,idir1,ipert1) 419 420 ! The first two dotprod corresponds to: <Psi(1)nkq|H(1)k+q,k|Psi(0)nk> and <Psi(0)nk|H(1)k,k+q|Psi(1)nkq> 421 ! They are calculated using wavefunctions <Psi(1)| that are orthogonal to the active space. 422 call dotprod_g(dotr,doti,istwf_k,npw1_k*nspinor,2,cwavef,gh1,mpi_enreg%me_g0,mpi_enreg%comm_spinorfft) 423 call dotprod_g(dot2r,dot2i,istwf_k,npw1_k*nspinor,2,gh,cwavef2,mpi_enreg%me_g0,mpi_enreg%comm_spinorfft) 424 425 ! This dotprod corresponds to : <Psi(1)nkq|H(0)k+q- E(0)nk|Psi(1)nkq> 426 ! It is calculated using wavefunctions that are orthogonal to the active space. 427 ! Should work for metals. (But adiabatic approximation is bad in this case...) 428 call dotprod_g(dot3r,dot3i,istwf_k,npw1_k*nspinor,2,cwavef,ghc,mpi_enreg%me_g0,mpi_enreg%comm_spinorfft) 429 430 eig2_diar= eig2_diar + dotr + dot2r + dot3r 431 eig2_diai= eig2_diai + doti + dot2i + dot3i 432 433 end if 434 435 ! Store the contribution 436 if(ieig2rf > 0) then 437 eig2nkq(1,iband+band_index,ikpt,idir1,ipert1,idir2,ipert2) = eig2_diar 438 eig2nkq(2,iband+band_index,ikpt,idir1,ipert1,idir2,ipert2) = eig2_diai 439 end if 440 441 if(present(eigbrd))then 442 if(smdelta >0) then !broadening 443 eigbrd(1,iband+band_index,ikpt,idir1,ipert1,idir2,ipert2) = eigbrd_r 444 eigbrd(2,iband+band_index,ikpt,idir1,ipert1,idir2,ipert2) = eigbrd_i 445 end if 446 end if 447 448 end do !idir2 449 end do !ipert2 450 end do !idir1 451 end do !ipert1 452 453 end if ! Selection of processor 454 455 end do !iband 456 457 ABI_DEALLOCATE(cwavef) 458 ABI_DEALLOCATE(cwavef2) 459 ABI_DEALLOCATE(gh) 460 ABI_DEALLOCATE(gh1) 461 ABI_DEALLOCATE(ghc) 462 band2tot_index = band2tot_index + 2*mband**2 463 bandtot_index = bandtot_index + mband 464 465 if(present(eigenq_fine))then 466 ABI_DEALLOCATE(kpt_fine_sub) ! Deallocate the variable 467 ABI_DEALLOCATE(wgt_sub) 468 end if 469 470 end do !ikpt 471 band_index = band_index + mband 472 end do !isppol 473 474 !Accumulate eig2nkq and/or eigbrd 475 if(xmpi_paral==1) then 476 if(ieig2rf == 1 .or. ieig2rf == 2) then 477 call xmpi_sum(eig2nkq,spaceworld,ierr) 478 if (dtset%dfpt_sciss > tol6 ) then 479 call xmpi_sum(eigen0,spaceworld,ierr) 480 call xmpi_sum(eigenq,spaceworld,ierr) 481 end if 482 end if 483 if(present(eigbrd) .and. (ieig2rf == 1 .or. ieig2rf == 2))then 484 if(smdelta >0) then 485 call xmpi_sum(eigbrd,spaceworld,ierr) 486 end if 487 end if 488 ABI_DEALLOCATE(nband_rbz) 489 ABI_DEALLOCATE(mpi_enreg%proc_distrb) 490 ABI_DEALLOCATE(mpi_enreg%my_kpttab) 491 end if 492 493 if(ieig2rf==1 .or. ieig2rf==2 ) then 494 write(ab_out,'(a)')' Components of second-order derivatives of the electronic energy, EIGR2D.' 495 write(ab_out,'(a)')' For automatic tests, printing the matrix for the first k-point, first band, first atom.' 496 do idir1=1,3 497 do idir2=1,3 498 ar=eig2nkq(1,1,1,idir1,1,idir2,1) ; if(abs(ar)<tol10)ar=zero 499 ai=eig2nkq(2,1,1,idir1,1,idir2,1) ; if(abs(ai)<tol10)ai=zero 500 write (ab_out,'(4i4,2es20.10)') idir1,1,idir2,1,ar,ai 501 end do ! idir2 502 end do ! idir1 503 end if 504 505 if(present(eigbrd))then 506 if(smdelta >0) then !broadening 507 write(ab_out,'(a)')' ' 508 write(ab_out,'(a)')' Components of second-order derivatives of the electronic energy, EIGI2D.' 509 write(ab_out,'(a)')' For automatic tests, printing the matrix for the first k-point, first band, first atom.' 510 do idir1=1,3 511 do idir2=1,3 512 ar=eigbrd(1,1,1,idir1,1,idir2,1) ; if(abs(ar)<tol10)ar=zero 513 ai=eigbrd(2,1,1,idir1,1,idir2,1) ; if(abs(ai)<tol10)ai=zero 514 write (ab_out,'(4i4,2es20.10)') idir1,1,idir2,1,ar,ai 515 end do 516 end do !nband 517 end if 518 end if 519 520 if(allocated(smdfun)) then 521 ABI_DEALLOCATE(smdfun) 522 end if 523 ABI_DEALLOCATE(icg2_rbz) 524 if(present(eigenq_fine))then 525 ABI_DEALLOCATE(center) 526 end if 527 if (dtset%dfpt_sciss > tol6 ) then 528 ABI_DEALLOCATE(eigen0tmp) 529 ABI_DEALLOCATE(eigenqtmp) 530 end if 531 532 call timab(148,2,tsec) 533 534 !DEBUG 535 !write(std_out,*)' eig2stern: exit' 536 !ENDDEBUG 537 538 end subroutine eig2stern