TABLE OF CONTENTS
ABINIT/eig2tot [ Functions ]
NAME
eig2tot
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 (SP,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.
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).
mband = maximum number of bands.
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.
nkpt_rbz = number of k-points for each perturbation.
nsppol = 1 for unpolarized, 2 for spin-polarized.
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 = header of the GS WF file of the corse k-point grid.
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
respfn
CHILDREN
crystal_free,crystal_init,ddb_hdr_free,ddb_hdr_init,ddb_hdr_open_write distrb2,ebands_free,ebands_init,eigr2d_free,eigr2d_init,eigr2d_ncwrite fan_free,fan_init,fan_ncwrite,gkk_free,gkk_init,gkk_ncwrite,kptfine_av outbsd,smeared_delta,timab,xmpi_sum
SOURCE
80 #if defined HAVE_CONFIG_H 81 #include "config.h" 82 #endif 83 84 #include "abi_common.h" 85 86 subroutine eig2tot(dtfil,xred,psps,pawtab,natom,bdeigrf,clflg,dim_eig2nkq,eigen0,eigenq,eigen1,eig2nkq,& 87 & elph2_imagden,esmear,ieig2rf,mband,mpert,npert,mpi_enreg,doccde,& 88 & nkpt_rbz,nsppol,smdelta,rprimd,dtset,occ_rbz,hdr0,eigbrd,eigenq_fine,hdr_fine) 89 90 91 use defs_basis 92 use defs_datatypes 93 use defs_abitypes 94 use m_profiling_abi 95 use m_xmpi 96 use m_nctk 97 use m_errors 98 #ifdef HAVE_NETCDF 99 use netcdf 100 #endif 101 use m_ebands 102 103 use m_fstrings, only : strcat 104 use m_eig2d, only : eigr2d_init,eigr2d_t,eigr2d_ncwrite,eigr2d_free,fan_t,& 105 & fan_init,fan_ncwrite,fan_free, gkk_t, gkk_init, & 106 & gkk_ncwrite,gkk_free 107 use m_crystal, only : crystal_init, crystal_free, crystal_t 108 use m_crystal_io, only : crystal_ncwrite 109 use m_pawtab, only : pawtab_type 110 use m_ddb, only : DDB_VERSION 111 use m_ddb_hdr, only : ddb_hdr_type, ddb_hdr_init, ddb_hdr_free, ddb_hdr_open_write 112 113 !This section has been created automatically by the script Abilint (TD). 114 !Do not modify the following lines by hand. 115 #undef ABI_FUNC 116 #define ABI_FUNC 'eig2tot' 117 use interfaces_18_timing 118 use interfaces_32_util 119 use interfaces_51_manage_mpi 120 use interfaces_72_response, except_this_one => eig2tot 121 !End of the abilint section 122 123 implicit none 124 125 !Arguments ------------------------------------ 126 !scalars 127 integer,intent(in) :: bdeigrf,dim_eig2nkq,ieig2rf,mband,mpert,natom,nkpt_rbz 128 integer,intent(in) :: npert,nsppol,smdelta 129 real(dp),intent(in) :: elph2_imagden,esmear 130 type(MPI_type),intent(inout) :: mpi_enreg 131 type(datafiles_type), intent(in) :: dtfil 132 type(pseudopotential_type), intent(inout) :: psps 133 !arrays 134 type(dataset_type), intent(in) :: dtset 135 integer,intent(in) :: clflg(3,mpert) 136 real(dp),intent(in) :: doccde(dtset%mband*dtset%nkpt*dtset%nsppol) 137 real(dp),intent(in) :: eigen0(nkpt_rbz*mband*nsppol) 138 real(dp),intent(in) :: eigen1(nkpt_rbz*2*nsppol*mband**2,3,mpert) 139 real(dp),intent(in) :: eigenq(nkpt_rbz*mband*nsppol) 140 real(dp),intent(in) :: occ_rbz(mband*nkpt_rbz*nsppol) 141 real(dp),intent(inout) :: eig2nkq(2,mband*nsppol,nkpt_rbz,3,npert,3,npert*dim_eig2nkq) 142 real(dp),intent(in) :: rprimd(3,3),xred(3,natom) 143 real(dp),intent(inout),optional :: eigbrd(2,mband*nsppol,nkpt_rbz,3,npert,3,npert) 144 real(dp),intent(in),pointer,optional :: eigenq_fine(:,:,:) 145 type(pawtab_type), intent(inout) :: pawtab(psps%ntypat*psps%usepaw) 146 type(hdr_type),intent(in) :: hdr0 147 type(hdr_type),intent(in),optional :: hdr_fine 148 149 !Local variables------------------------------- 150 !tolerance for non degenerated levels 151 !scalars 152 integer :: band2tot_index,band_index,bantot,bandtot_index,iband,idir1,idir2 153 integer :: ikpt,ipert1,ipert2,isppol,jband,nkpt_sub,ikpt2,unitout,ncid 154 !integer :: ipw 155 character(len=fnlen) :: dscrpt,fname 156 integer :: master,me,spaceworld,ierr 157 ! real(dp),parameter :: etol=1.0d-6 158 real(dp),parameter :: etol=1.0d-7 159 !real(dp),parameter :: etol=zero 160 real(dp) :: ar,ai,deltae,den,eig1_i1,eig1_i2,eigen_corr 161 real(dp) :: eig1_r1,eig1_r2,eig2_diai,den_av 162 real(dp) :: eig2_diar,eigbrd_i,eigbrd_r,wgt_int 163 character(len=500) :: message 164 logical :: remove_inv,test_do_band 165 type(crystal_t) :: Crystal 166 type(ebands_t) :: Bands 167 type(eigr2d_t) :: eigr2d,eigi2d 168 type(fan_t) :: fan2d 169 type(gkk_t) :: gkk2d 170 type(ddb_hdr_type) :: ddb_hdr 171 !arrays 172 integer, allocatable :: nband_rbz(:) 173 integer,pointer :: kpt_fine_sub(:) 174 real(dp) :: tsec(2) 175 real(dp),allocatable :: center(:) 176 real(dp) :: eigen(mband*nsppol),eigen_prime(mband*nsppol) 177 real(dp),allocatable :: fan(:,:,:,:,:,:,:) 178 real(dp),allocatable :: gkk(:,:,:,:,:) 179 real(dp),allocatable :: eig2nkq_tmp(:,:,:,:,:,:,:) 180 real(dp),allocatable :: smdfun(:,:) 181 real(dp),pointer :: wgt_sub(:) 182 183 ! ********************************************************************* 184 185 !Init parallelism 186 master =0 187 spaceworld=mpi_enreg%comm_cell 188 me=mpi_enreg%me_kpt 189 !DEBUG 190 !write(std_out,*)' eig2tot : enter ' 191 !write(std_out,*)' mband=',mband 192 !write(std_out,*)' nsppol=',nsppol 193 !write(std_out,*)' nkpt_rbz=',nkpt_rbz 194 !write(std_out,*)' npert=',npert 195 !ENDDEBUG 196 197 !Init interpolation method 198 if(present(eigenq_fine))then 199 ABI_ALLOCATE(center,(3)) 200 end if 201 202 call timab(148,1,tsec) 203 204 if(nsppol==2)then 205 message = 'nsppol=2 is still under development. Be careful when using it ...' 206 MSG_COMMENT(message) 207 end if 208 209 band2tot_index =0 210 bandtot_index=0 211 band_index=0 212 213 if(xmpi_paral==1) then 214 ABI_ALLOCATE(mpi_enreg%proc_distrb,(nkpt_rbz,mband,nsppol)) 215 ABI_ALLOCATE(nband_rbz,(nkpt_rbz*nsppol)) 216 if (allocated(mpi_enreg%my_kpttab)) then 217 ABI_DEALLOCATE(mpi_enreg%my_kpttab) 218 end if 219 ABI_ALLOCATE(mpi_enreg%my_kpttab,(nkpt_rbz)) 220 ! Assume the number of bands is the same for all k points. 221 nband_rbz(:)=mband 222 call distrb2(mband,nband_rbz,nkpt_rbz,mpi_enreg%nproc_cell,nsppol,mpi_enreg) 223 end if 224 225 if(ieig2rf == 4 ) then 226 ABI_STAT_ALLOCATE(fan,(2*mband*nsppol,dtset%nkpt,3,natom,3,natom*dim_eig2nkq,mband), ierr) 227 ABI_CHECK(ierr==0, "out-of-memory in fan") 228 fan(:,:,:,:,:,:,:) = zero 229 ABI_STAT_ALLOCATE(eig2nkq_tmp,(2,mband*nsppol,dtset%nkpt,3,natom,3,natom*dim_eig2nkq), ierr) 230 ABI_CHECK(ierr==0, "out-of-memory in eig2nkq_tmp") 231 eig2nkq_tmp(:,:,:,:,:,:,:) = zero 232 ! This is not efficient because double the memory. Alternative: use buffer and 233 ! print part by part. 234 eig2nkq_tmp = eig2nkq 235 if(present(eigbrd))then 236 eigbrd(:,:,:,:,:,:,:)=zero 237 end if 238 eigen_corr = 0 239 end if 240 241 if(ieig2rf == 5 ) then 242 ABI_STAT_ALLOCATE(gkk,(2*mband*nsppol,dtset%nkpt,3,natom,mband), ierr) 243 ABI_CHECK(ierr==0, "out-of-memory in gkk") 244 gkk(:,:,:,:,:) = zero 245 ABI_STAT_ALLOCATE(eig2nkq_tmp,(2,mband*nsppol,dtset%nkpt,3,natom,3,natom*dim_eig2nkq), ierr) 246 ABI_CHECK(ierr==0, "out-of-memory in eig2nkq_tmp") 247 eig2nkq_tmp(:,:,:,:,:,:,:) = zero 248 ! This is not efficient because double the memory. Alternative: use buffer and 249 ! print part by part. 250 eig2nkq_tmp = eig2nkq 251 if(present(eigbrd))then 252 eigbrd(:,:,:,:,:,:,:)=zero 253 end if 254 eigen_corr = 0 255 end if 256 257 do isppol=1,nsppol 258 do ikpt =1,nkpt_rbz 259 260 if(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,mband,isppol,me)) then 261 band2tot_index = band2tot_index + 2*mband**2 262 bandtot_index = bandtot_index + mband 263 cycle 264 end if 265 266 if(present(eigenq_fine))then 267 write(std_out,*) 'Start of the energy denominator interpolation method.' 268 nkpt_sub = 0 269 ! center is the k+q point around which we will average the kpt_fine 270 center = hdr0%kptns(:,ikpt)+ dtset%qptn(:) 271 272 call kptfine_av(center,dtset%qptrlatt,hdr_fine%kptns,hdr_fine%nkpt,& 273 & kpt_fine_sub,nkpt_sub,wgt_sub) 274 write(std_out,'(a,3f8.4,a,i3)') 'Number of k-points of the fine grid & 275 & around the k+Q point ',center,' is:',nkpt_sub 276 write(std_out,'(a,f10.5)') 'The sum of the weights of the k-points is: ',SUM(wgt_sub) 277 end if 278 279 if(smdelta >0) then !broadening 280 if(.not.allocated(smdfun)) then 281 ABI_ALLOCATE(smdfun,(mband,mband)) 282 end if 283 smdfun(:,:) = zero 284 do iband=1,mband 285 eigen(iband) = eigen0(iband+bandtot_index) 286 eigen_prime(iband) =eigenq(iband+bandtot_index) 287 end do 288 if(esmear>tol6) then 289 call smeared_delta(eigen,eigen_prime,esmear,mband,smdelta,smdfun) 290 end if 291 end if 292 293 ipert1=1 ! Suppose all perturbations lead to the same number of planewaves 294 295 do iband=1,bdeigrf 296 297 ! If the k point and band belong to me, compute the contribution 298 test_do_band=.true. 299 if(mpi_enreg%proc_distrb(ikpt,iband,isppol)/=me)test_do_band=.false. 300 301 if(test_do_band)then 302 ! ------------------------------------------------------------------------------------------------------! 303 ! ------- ieig2rf ==3 : Non dynamic traditional AHC theory with Sternheimer (computed in eig2stern.F90)-! 304 ! ------------------------------------------------------------------------------------------------------! 305 ! Note that ieig2rf==4 and ieig2rf==5 also goes into that part only for later printing of the ZPR in the ouput of abinit 306 ! later in the code 307 if(ieig2rf==3 .or. ieig2rf==4 .or. ieig2rf==5) then 308 do ipert1=1,npert 309 do idir1=1,3 310 if(clflg(idir1,ipert1)==0) cycle 311 do ipert2=1,npert 312 do idir2=1,3 313 if(clflg(idir2,ipert2)==0)cycle 314 eig2_diar = zero ; eig2_diai = zero ; eigbrd_r = zero ; eigbrd_i = zero 315 do jband=1,mband 316 eig1_r1 = eigen1(2*jband-1+(iband-1)*2*mband+band2tot_index,idir1,ipert1) 317 eig1_r2 = eigen1(2*jband-1+(iband-1)*2*mband+band2tot_index,idir2,ipert2) 318 eig1_i1 = eigen1(2*jband+(iband-1)*2*mband+band2tot_index,idir1,ipert1) 319 eig1_i2 = - eigen1(2*jband+(iband-1)*2*mband+band2tot_index,idir2,ipert2) !the negative sign is from the CC 320 ! If no interpolation, fallback on to the previous 321 ! implementation 322 if(.not. present(eigenq_fine))then 323 deltae=eigenq(jband+bandtot_index)-eigen0(iband+bandtot_index) 324 end if 325 ar=eig1_r1*eig1_r2-eig1_i1*eig1_i2 326 ai=eig1_r1*eig1_i2+eig1_i1*eig1_r2 327 328 ! Sum over all active space to retrieve the diagonal gauge 329 ! if(abs(deltae)>etol) then ! This is commented because 330 ! there is no problem with divergencies with elph2_imag != 0 331 if( present(eigenq_fine))then 332 den_av = zero 333 wgt_int = zero 334 do ikpt2=1,nkpt_sub 335 deltae=eigenq_fine(jband,kpt_fine_sub(ikpt2),1)& 336 & -eigen0(iband+bandtot_index) 337 den_av = den_av-(wgt_sub(ikpt2)*deltae)/(deltae**2+elph2_imagden**2) 338 wgt_int = wgt_int+wgt_sub(ikpt2) 339 end do 340 den = den_av/wgt_int 341 else 342 if(abs(elph2_imagden) < etol) then 343 if(abs(deltae)>etol) then 344 den=-one/(deltae**2+elph2_imagden**2) 345 else 346 den= zero 347 end if 348 else 349 den=-one/(deltae**2+elph2_imagden**2) 350 end if 351 end if 352 353 if( present(eigenq_fine))then 354 eig2_diar=eig2_diar+ar*den 355 eig2_diai=eig2_diai+ai*den 356 else 357 eig2_diar=eig2_diar+ar*deltae*den 358 eig2_diai=eig2_diai+ai*deltae*den 359 end if 360 361 if(present(eigbrd))then 362 if(smdelta >0) then !broadening 363 eigbrd_r = eigbrd_r + ar*smdfun(iband,jband) 364 eigbrd_i = eigbrd_i + ai*smdfun(iband,jband) 365 end if 366 end if 367 end do !jband 368 369 ! Store the contribution 370 eig2nkq(1,iband+band_index,ikpt,idir1,ipert1,idir2,ipert2) = & 371 & eig2nkq(1,iband+band_index,ikpt,idir1,ipert1,idir2,ipert2) + eig2_diar 372 eig2nkq(2,iband+band_index,ikpt,idir1,ipert1,idir2,ipert2) = & 373 & eig2nkq(2,iband+band_index,ikpt,idir1,ipert1,idir2,ipert2) + eig2_diai 374 375 if(present(eigbrd))then 376 if(smdelta >0) then !broadening 377 eigbrd(1,iband+band_index,ikpt,idir1,ipert1,idir2,ipert2) = eigbrd_r 378 eigbrd(2,iband+band_index,ikpt,idir1,ipert1,idir2,ipert2) = eigbrd_i 379 end if 380 end if 381 382 end do !idir2 383 end do !ipert2 384 end do !idir1 385 end do !ipert1 386 end if !ieig2rf 3 387 388 ! -------------------------------------------------------------------------------------------! 389 ! ------- ieig2rf ==4 Dynamic AHC using second quantization and Sternheimer from eig2stern -! 390 ! -------------------------------------------------------------------------------------------! 391 if(ieig2rf ==4 ) then 392 do ipert1=1,npert 393 do idir1=1,3 394 if(clflg(idir1,ipert1)==0) cycle 395 do ipert2=1,npert 396 do idir2=1,3 397 if(clflg(idir2,ipert2)==0)cycle 398 do jband=1,mband 399 eig1_r1 = eigen1(2*jband-1+(iband-1)*2*mband+band2tot_index,idir1,ipert1) 400 eig1_r2 = eigen1(2*jband-1+(iband-1)*2*mband+band2tot_index,idir2,ipert2) 401 eig1_i1 = eigen1(2*jband+(iband-1)*2*mband+band2tot_index,idir1,ipert1) 402 eig1_i2 = - eigen1(2*jband+(iband-1)*2*mband+band2tot_index,idir2,ipert2) !the negative sign is from the CC 403 ar=eig1_r1*eig1_r2-eig1_i1*eig1_i2 404 ai=eig1_r1*eig1_i2+eig1_i1*eig1_r2 405 ! Store the contribution 406 fan(2*iband-1+2*band_index,ikpt,idir1,ipert1,idir2,ipert2,jband) = & 407 & fan(2*iband-1+2*band_index,ikpt,idir1,ipert1,idir2,ipert2,jband) + ar 408 fan(2*iband+2*band_index,ikpt,idir1,ipert1,idir2,ipert2,jband) = & 409 & fan(2*iband+2*band_index,ikpt,idir1,ipert1,idir2,ipert2,jband) + ai 410 end do !jband 411 end do !idir2 412 end do !ipert2 413 end do !idir1 414 end do !ipert1 415 end if !ieig2rf 4 416 ! --------------------------------------------------------------------------------! 417 ! ------- ieig2rf ==5 Dynamic AHC with Sternheimer from eig2stern but print GKK -! 418 ! --------------------------------------------------------------------------------! 419 if(ieig2rf ==5 ) then 420 do ipert1=1,npert 421 do idir1=1,3 422 if(clflg(idir1,ipert1)==0) cycle 423 do jband=1,mband 424 eig1_r1 = eigen1(2*jband-1+(iband-1)*2*mband+band2tot_index,idir1,ipert1) 425 eig1_i1 = eigen1(2*jband+(iband-1)*2*mband+band2tot_index,idir1,ipert1) 426 ! Store the contribution 427 gkk(2*iband-1+2*band_index,ikpt,idir1,ipert1,jband) = & 428 & gkk(2*iband-1+2*band_index,ikpt,idir1,ipert1,jband) + eig1_r1 429 gkk(2*iband+2*band_index,ikpt,idir1,ipert1,jband) = & 430 & gkk(2*iband+2*band_index,ikpt,idir1,ipert1,jband) + eig1_i1 431 end do !jband 432 end do !idir1 433 end do !ipert1 434 end if !ieig2rf 5 435 end if ! Selection of processor 436 end do !iband 437 438 band2tot_index = band2tot_index + 2*mband**2 439 bandtot_index = bandtot_index + mband 440 441 if(present(eigenq_fine))then 442 ABI_DEALLOCATE(kpt_fine_sub) ! Deallocate the variable 443 ABI_DEALLOCATE(wgt_sub) 444 end if 445 446 end do !ikpt 447 band_index = band_index + mband 448 end do !isppol 449 450 !Accumulate eig2nkq and/or eigbrd 451 if(xmpi_paral==1) then 452 if(ieig2rf == 3) then 453 call xmpi_sum(eig2nkq,spaceworld,ierr) 454 end if 455 if(ieig2rf == 4) then 456 call xmpi_sum(eig2nkq,spaceworld,ierr) 457 call xmpi_sum(eig2nkq_tmp,spaceworld,ierr) 458 call xmpi_sum(fan,spaceworld,ierr) 459 end if 460 if(ieig2rf == 5) then 461 call xmpi_sum(eig2nkq,spaceworld,ierr) 462 call xmpi_sum(eig2nkq_tmp,spaceworld,ierr) 463 call xmpi_sum(gkk,spaceworld,ierr) 464 end if 465 if(present(eigbrd) .and. (ieig2rf == 3 .or. ieig2rf == 4 .or. ieig2rf == 5))then 466 if(smdelta >0) then 467 call xmpi_sum(eigbrd,spaceworld,ierr) 468 end if 469 end if 470 ABI_DEALLOCATE(nband_rbz) 471 ABI_DEALLOCATE(mpi_enreg%proc_distrb) 472 ABI_DEALLOCATE(mpi_enreg%my_kpttab) 473 end if 474 475 if(ieig2rf > 2) then 476 write(ab_out,'(a)')' Components of second-order derivatives of the electronic energy, EIGR2D.' 477 write(ab_out,'(a)')' For automatic tests, printing the matrix for the first k-point, first band, first atom.' 478 band_index = 0 479 do isppol=1,dtset%nsppol 480 do idir1=1,3 481 do idir2=1,3 482 ar=eig2nkq(1,1+band_index,1,idir1,1,idir2,1) ; if(abs(ar)<tol10)ar=zero 483 ai=eig2nkq(2,1+band_index,1,idir1,1,idir2,1) ; if(abs(ai)<tol10)ai=zero 484 write (ab_out,'(4i4,2es20.10)') idir1,1,idir2,1,ar,ai 485 end do ! idir2 486 end do ! idir1 487 band_index = band_index + mband 488 write(ab_out,'(a)')' ' 489 end do 490 end if 491 492 if(present(eigbrd))then 493 if(smdelta >0) then !broadening 494 write(ab_out,'(a)')' Components of second-order derivatives of the electronic energy, EIGI2D.' 495 write(ab_out,'(a)')' For automatic tests, printing the matrix for the first k-point, first band, first atom.' 496 band_index = 0 497 do isppol=1,dtset%nsppol 498 do idir1=1,3 499 do idir2=1,3 500 ar=eigbrd(1,1+band_index,1,idir1,1,idir2,1) ; if(abs(ar)<tol10)ar=zero 501 ai=eigbrd(2,1+band_index,1,idir1,1,idir2,1) ; if(abs(ai)<tol10)ai=zero 502 write (ab_out,'(4i4,2es20.10)') idir1,1,idir2,1,ar,ai 503 end do 504 end do 505 band_index = band_index + mband 506 write(ab_out,'(a)')' ' 507 end do 508 end if 509 end if 510 511 if(allocated(smdfun)) then 512 ABI_DEALLOCATE(smdfun) 513 end if 514 if(present(eigenq_fine))then 515 ABI_DEALLOCATE(center) 516 end if 517 518 master=0 519 if (me==master) then 520 ! print _EIGR2D file for this perturbation in the case of ieig2rf 3 or 4 or 5 521 if (ieig2rf == 3 .or. ieig2rf == 4 .or. ieig2rf == 5) then 522 523 dscrpt=' Note : temporary (transfer) database ' 524 unitout = dtfil%unddb 525 526 call ddb_hdr_init(ddb_hdr,dtset,psps,pawtab,DDB_VERSION,dscrpt,& 527 & 1,xred=xred,occ=occ_rbz) 528 529 call ddb_hdr_open_write(ddb_hdr, dtfil%fnameabo_eigr2d, unitout) 530 531 call ddb_hdr_free(ddb_hdr) 532 533 end if 534 if(ieig2rf == 3 ) then 535 call outbsd(bdeigrf,dtset,eig2nkq,dtset%natom,nkpt_rbz,unitout) 536 end if 537 if(ieig2rf == 4 .or. ieig2rf == 5 ) then 538 call outbsd(bdeigrf,dtset,eig2nkq_tmp,dtset%natom,nkpt_rbz,unitout) 539 end if 540 ! Output of the EIGR2D.nc file. 541 fname = strcat(dtfil%filnam_ds(4),"_EIGR2D.nc") 542 ! Crystalline structure. 543 remove_inv=.false. 544 if(dtset%nspden==4 .and. dtset%usedmft==1) remove_inv=.true. 545 call crystal_init(dtset%amu_orig(:,1),Crystal,dtset%spgroup,dtset%natom,dtset%npsp,psps%ntypat, & 546 & dtset%nsym,rprimd,dtset%typat,xred,dtset%ziontypat,dtset%znucl,1,& 547 & dtset%nspden==2.and.dtset%nsppol==1,remove_inv,hdr0%title,& 548 & dtset%symrel,dtset%tnons,dtset%symafm) 549 ! Electronic band energies. 550 bantot= dtset%mband*dtset%nkpt*dtset%nsppol 551 call ebands_init(bantot,Bands,dtset%nelect,doccde,eigen0,hdr0%istwfk,hdr0%kptns,& 552 & hdr0%nband, hdr0%nkpt,hdr0%npwarr,hdr0%nsppol,hdr0%nspinor,& 553 & hdr0%tphysel,hdr0%tsmear,hdr0%occopt,hdr0%occ,hdr0%wtk,& 554 & hdr0%charge, hdr0%kptopt, hdr0%kptrlatt_orig, hdr0%nshiftk_orig, hdr0%shiftk_orig, & 555 & hdr0%kptrlatt, hdr0%nshiftk, hdr0%shiftk) 556 557 ! Second order derivative EIGR2D (real and Im) 558 if(ieig2rf == 3 ) then 559 call eigr2d_init(eig2nkq,eigr2d,dtset%mband,hdr0%nsppol,nkpt_rbz,dtset%natom) 560 end if 561 if(ieig2rf == 4 .or. ieig2rf == 5 ) then 562 call eigr2d_init(eig2nkq_tmp,eigr2d,dtset%mband,hdr0%nsppol,nkpt_rbz,dtset%natom) 563 end if 564 #ifdef HAVE_NETCDF 565 NCF_CHECK_MSG(nctk_open_create(ncid, fname, xmpi_comm_self), "Creating EIGR2D file") 566 NCF_CHECK(crystal_ncwrite(Crystal,ncid)) 567 NCF_CHECK(ebands_ncwrite(Bands, ncid)) 568 call eigr2d_ncwrite(eigr2d,dtset%qptn(:),dtset%wtq,ncid) 569 NCF_CHECK(nf90_close(ncid)) 570 #else 571 ABI_UNUSED(ncid) 572 #endif 573 574 ! print _FAN file for this perturbation. Note that the Fan file will only be produced if 575 ! abinit is compiled with netcdf. 576 if(ieig2rf == 4 ) then 577 ! Output of the Fan.nc file. 578 #ifdef HAVE_NETCDF 579 fname = strcat(dtfil%filnam_ds(4),"_FAN.nc") 580 call fan_init(fan,fan2d,dtset%mband,hdr0%nsppol,nkpt_rbz,dtset%natom) 581 NCF_CHECK_MSG(nctk_open_create(ncid, fname, xmpi_comm_self), "Creating FAN file") 582 NCF_CHECK(crystal_ncwrite(Crystal, ncid)) 583 NCF_CHECK(ebands_ncwrite(Bands, ncid)) 584 call fan_ncwrite(fan2d,dtset%qptn(:),dtset%wtq, ncid) 585 NCF_CHECK(nf90_close(ncid)) 586 #else 587 MSG_ERROR("Dynamical calculation with ieig2rf 4 only work with NETCDF support.") 588 ABI_UNUSED(ncid) 589 #endif 590 ABI_DEALLOCATE(fan) 591 ABI_DEALLOCATE(eig2nkq_tmp) 592 end if 593 ! print _GKK.nc file for this perturbation. Note that the GKK file will only be produced if 594 ! abinit is compiled with netcdf. 595 if(ieig2rf == 5 ) then 596 ! Output of the GKK.nc file. 597 #ifdef HAVE_NETCDF 598 fname = strcat(dtfil%filnam_ds(4),"_GKK.nc") 599 call gkk_init(gkk,gkk2d,dtset%mband,hdr0%nsppol,nkpt_rbz,dtset%natom,3) 600 NCF_CHECK_MSG(nctk_open_create(ncid, fname, xmpi_comm_self), "Creating GKK file") 601 NCF_CHECK(crystal_ncwrite(Crystal, ncid)) 602 NCF_CHECK(ebands_ncwrite(Bands, ncid)) 603 call gkk_ncwrite(gkk2d,dtset%qptn(:),dtset%wtq, ncid) 604 NCF_CHECK(nf90_close(ncid)) 605 #else 606 MSG_ERROR("Dynamical calculation with ieig2rf 5 only work with NETCDF support.") 607 ABI_UNUSED(ncid) 608 #endif 609 ABI_DEALLOCATE(gkk) 610 ABI_DEALLOCATE(eig2nkq_tmp) 611 end if 612 ! print _EIGI2D file for this perturbation 613 if (ieig2rf /= 5 ) then 614 if(smdelta>0) then 615 unitout = dtfil%unddb 616 dscrpt=' Note : temporary (transfer) database ' 617 618 call ddb_hdr_init(ddb_hdr,dtset,psps,pawtab,DDB_VERSION,dscrpt,& 619 & 1,xred=xred,occ=occ_rbz) 620 621 call ddb_hdr_open_write(ddb_hdr, dtfil%fnameabo_eigi2d, unitout) 622 623 call ddb_hdr_free(ddb_hdr) 624 625 call outbsd(bdeigrf,dtset,eigbrd,dtset%natom,nkpt_rbz,unitout) 626 627 ! Output of the EIGI2D.nc file. 628 fname = strcat(dtfil%filnam_ds(4),"_EIGI2D.nc") 629 ! Broadening EIGI2D (real and Im) 630 call eigr2d_init(eigbrd,eigi2d,dtset%mband,hdr0%nsppol,nkpt_rbz,dtset%natom) 631 #ifdef HAVE_NETCDF 632 NCF_CHECK_MSG(nctk_open_create(ncid, fname, xmpi_comm_self), "Creating EIGI2D file") 633 NCF_CHECK(crystal_ncwrite(Crystal, ncid)) 634 NCF_CHECK(ebands_ncwrite(Bands, ncid)) 635 call eigr2d_ncwrite(eigi2d,dtset%qptn(:),dtset%wtq,ncid) 636 NCF_CHECK(nf90_close(ncid)) 637 #else 638 ABI_UNUSED(ncid) 639 #endif 640 end if !smdelta 641 end if 642 end if 643 644 if (allocated(fan)) then 645 ABI_DEALLOCATE(fan) 646 end if 647 if (allocated(eig2nkq_tmp)) then 648 ABI_DEALLOCATE(eig2nkq_tmp) 649 end if 650 if (allocated(gkk)) then 651 ABI_DEALLOCATE(gkk) 652 end if 653 654 call crystal_free(Crystal) 655 call ebands_free(Bands) 656 call eigr2d_free(eigr2d) 657 call eigr2d_free(eigi2d) 658 call fan_free(fan2d) 659 call gkk_free(gkk2d) 660 661 662 call timab(148,2,tsec) 663 !DEBUG 664 !write(std_out,*)' eig2tot: exit' 665 !ENDDEBUG 666 667 end subroutine eig2tot