TABLE OF CONTENTS
ABINIT/optic [ Programs ]
NAME
optic
FUNCTION
Driver routine to call linopt and nlinopt, which calculate the linear and non-linear optical responses in the RPA.
COPYRIGHT
Copyright (C) 2002-2024 ABINIT group (SSharma,MVer,VRecoules,YG,NAP) 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
(main routine)
OUTPUT
(main routine)
NOTES
domega=frequency range eigen11(2*mband*mband*nkpt_rbz*nsppol)=first-order eigenvalues (hartree) in reciprocal direction 100 eigen12(2*mband*mband*nkpt_rbz*nsppol)=first-order eigenvalues (hartree) in reciprocal direction 010 eigen13(2*mband*mband*nkpt_rbz*nsppol)=first-order eigenvalues (hartree) in reciprocal direction 001 nomega=number of frequency for conductivity computation mband=maximum number of bands. occopt==option for occupancies broadening=smearing width (or temperature) in Hartree maxomega=frequency windows for computations of sigma
SOURCE
35 #if defined HAVE_CONFIG_H 36 #include "config.h" 37 #endif 38 39 #include "abi_common.h" 40 41 program optic 42 43 use defs_basis 44 use m_errors 45 use m_xmpi 46 use m_xomp 47 use m_abicore 48 use m_optic_tools 49 use m_wfk 50 use m_nctk 51 use m_hdr 52 use m_ebands 53 use m_eprenorms 54 use m_crystal 55 use m_argparse 56 #ifdef HAVE_NETCDF 57 use netcdf 58 #endif 59 60 use m_build_info, only : abinit_version 61 use defs_datatypes, only : ebands_t 62 use m_specialmsg, only : specialmsg_getcount, herald 63 use m_time , only : asctime, timein 64 use m_symtk, only : mati3inv, matr3inv 65 use m_geometry, only : metric 66 use m_io_tools, only : flush_unit, open_file, file_exists, get_unit 67 use m_numeric_tools, only : c2r 68 use m_fstrings, only : int2char4, itoa, sjoin, strcat, endswith, basename 69 70 implicit none 71 72 !Arguments ----------------------------------- 73 74 !Local variables------------------------------- 75 integer,parameter :: formeig0 = 0, formeig1 = 1, master = 0 76 integer :: fform,finunt,ep_ntemp,itemp,i1,i2 77 integer :: bantot,bdtot0_index,bdtot_index 78 integer :: ierr,ii,jj,ikpt 79 integer :: isppol,mband,nomega,nband1 80 integer :: nkpt,nsppol 81 integer :: nks_per_proc,work_size,lin1,lin2,nlin1,nlin2,nlin3 82 integer :: linel1,linel2,linel3,nonlin1,nonlin2,nonlin3 83 integer :: iomode0,comm,nproc,my_rank, optic_ncid 84 #ifdef HAVE_NETCDF 85 integer :: ncid, varid, ncerr 86 #endif 87 integer :: num_lin_comp=1,num_nonlin_comp=0,num_linel_comp=0,num_nonlin2_comp=0 88 integer :: autoparal=0,max_ncpus=0 89 integer :: nonlin_comp(27) = 0, linel_comp(27) = 0, nonlin2_comp(27) = 0 90 integer :: lin_comp(9) = [11, 22 ,33, 12, 13, 21, 23, 31, 32] 91 integer :: prtlincompmatrixelements=0, nband_sum = -1 92 real(dp) :: domega, eff 93 real(dp) :: broadening,maxomega,scissor,tolerance 94 real(dp) :: tcpu,tcpui,twall,twalli 95 logical :: do_antiresonant, do_temperature 96 logical :: do_ep_renorm 97 logical,parameter :: remove_inv = .False. 98 type(hdr_type) :: hdr 99 type(ebands_t) :: ks_ebands, eph_ebands 100 type(crystal_t) :: cryst 101 type(eprenorms_t) :: Epren 102 type(wfk_t) :: wfk0 103 type(args_t) :: args 104 !arrays 105 integer :: iomode_ddk(3) 106 real(dp) :: tsec(2) 107 real(dp),allocatable :: wmesh(:) 108 real(dp),allocatable :: doccde(:), eig0tmp(:), eigen0(:) 109 real(dp),target,allocatable :: eigen11(:),eigen12(:),eigen13(:) 110 real(dp),allocatable :: eigtmp(:) 111 real(dp), ABI_CONTIGUOUS pointer :: outeig(:) 112 complex(dpc),allocatable :: pmat(:,:,:,:,:) 113 logical :: use_ncevk(0:3) 114 character(len=fnlen) :: filnam,wfkfile,ddkfile_1,ddkfile_2,ddkfile_3,filnam_out, epfile,fname 115 character(len=fnlen) :: infiles(0:3) 116 character(len=256) :: prefix,tmp_radix 117 character(len=10) :: s1,s2,s3,stemp 118 character(len=24) :: codename, start_datetime 119 character(len=500) :: msg 120 character(len=fnlen) :: ep_nc_fname 121 type(hdr_type) :: hdr_ddk(3) 122 type(wfk_t) :: wfks(0:3) 123 124 ! Input file 125 namelist /FILES/ ddkfile_1, ddkfile_2, ddkfile_3, wfkfile 126 namelist /PARAMETERS/ broadening, domega, maxomega, scissor, tolerance, do_antiresonant, do_temperature, & 127 autoparal, max_ncpus, prtlincompmatrixelements, nband_sum 128 namelist /COMPUTATIONS/ num_lin_comp, lin_comp, num_nonlin_comp, nonlin_comp, & 129 num_linel_comp, linel_comp, num_nonlin2_comp, nonlin2_comp 130 namelist /TEMPERATURE/ epfile 131 132 ! ********************************************************************************* 133 134 ! Change communicator for I/O (mandatory!) 135 call abi_io_redirect(new_io_comm=xmpi_world) 136 137 call xmpi_init() 138 comm = xmpi_world 139 nproc = xmpi_comm_size(comm); my_rank = xmpi_comm_rank(comm) 140 141 ! Parse command line arguments. 142 args = args_parser(); if (args%exit /= 0) goto 100 143 144 ! Initialize memory profiling if it is activated 145 ! if a full abimem.mocc report is desired, set the argument of abimem_init to "2" instead of "0" 146 ! note that abimem.mocc files can easily be multiple GB in size so don't use this option normally 147 #ifdef HAVE_MEM_PROFILING 148 call abimem_init(args%abimem_level, limit_mb=args%abimem_limit_mb) 149 #endif 150 151 call timein(tcpui,twall) 152 call timein(tcpui,twalli) 153 start_datetime = asctime() 154 155 if (my_rank == master) then 156 codename='OPTIC '//repeat(' ',18) 157 call herald(codename,abinit_version,std_out) 158 159 if (len_trim(args%input_path) == 0) then 160 ! Legacy Files file mode. 161 write(std_out, "(2a)")" DeprecationWarning: ",ch10 162 write(std_out, "(a)") " The files file has been deprecated in Abinit9 and will be removed in Abinit10." 163 write(std_out, "(a)")" Use the syntax `optic t01.abi` to run optic" 164 165 !Read data file name 166 write(std_out,'(a)')' Please, give the name of the data file ...' 167 read(5, '(a)')filnam 168 write(std_out,'(a,a,1x,a,a)')' The name of the data file is :',ch10,trim(filnam),ch10 169 write(std_out,'(a)')' Please, give the name of the output file ...' 170 read(5, '(a)')filnam_out 171 write(std_out,'(a,a,1x,a,a)')' The name of the output file is :',ch10,trim(filnam_out),ch10 172 write(std_out,'(a)')' Please, give the root name for the (non)linear optical data output file ...' 173 read(5, '(a)')prefix 174 write(std_out,'(a,a,1x,a)')' The root name of the output files is :',ch10,trim(prefix) 175 176 else 177 filnam = args%input_path 178 ! Get prefix from input file. Default values are provided 179 filnam_out = trim(filnam)//".abo" 180 prefix = trim(filnam) 181 182 ! If the basename has file extension e.g. run.abi, use what comes before the dot to build 183 ! filnam_out (e.g. run.abo) and the prefix for output files. 184 fname = basename(args%input_path) 185 i1 = index(fname, ".") 186 if (i1 > 1) then 187 i2 = index(args%input_path, ".", back=.True.) 188 filnam_out = args%input_path(:i2) // "abo" 189 prefix = args%input_path(:i2-1) 190 end if 191 end if 192 193 ! Read data file 194 if (open_file(filnam,msg,newunit=finunt,form='formatted') /= 0) then 195 ABI_ERROR(msg) 196 end if 197 198 ! Setup some default values: 199 broadening = 1e-3_dp ! Ha 200 domega = 1e-3_dp ! Ha 201 maxomega = 1.0_dp ! Ha 202 scissor = 0.0_dp ! no scissor by default 203 tolerance = 1e-3_dp ! Ha 204 prtlincompmatrixelements = 0 ! print the sum elements for external analysis 205 do_antiresonant = .TRUE. ! do use antiresonant approximation (only resonant transitions in the calculation) 206 do_temperature = .FALSE. 207 208 ! Read input file 209 read(finunt,nml=FILES) 210 read(finunt,nml=PARAMETERS) 211 read(finunt,nml=COMPUTATIONS) 212 if (do_temperature) read(finunt, nml=TEMPERATURE) 213 close(finunt) 214 ! Store filenames in array. 215 infiles = [wfkfile, ddkfile_1, ddkfile_2, ddkfile_3] 216 217 ! Validate input 218 if (num_nonlin_comp > 0 .and. all(nonlin_comp(1:num_nonlin_comp) == 0)) then 219 ABI_ERROR("nonlin_comp must be specified when num_nonlin_comp > 0") 220 end if 221 if (num_linel_comp > 0 .and. all(linel_comp(1:num_linel_comp) == 0)) then 222 ABI_ERROR("linel_comp must be specified when num_linel_comp > 0") 223 end if 224 if (num_nonlin2_comp > 0 .and. all(nonlin2_comp(1:num_nonlin2_comp) == 0)) then 225 ABI_ERROR("nonlin2_comp must be specified when num_nonlin2_comp > 0") 226 end if 227 228 ! Open GS wavefunction file 229 ! Note: Cannot use MPI-IO here because of prtwf=3. 230 ! If prtwf==3, the DDK file does not contain the wavefunctions but 231 ! this info is not reported in the header and the offsets in wfk_compute_offsets 232 ! are always computed assuming the presence of the cg 233 call nctk_fort_or_ncfile(wfkfile, iomode0, msg) 234 if (len_trim(msg) /= 0) ABI_ERROR(msg) 235 if (iomode0 == IO_MODE_MPI) iomode0 = IO_MODE_FORTRAN 236 call wfk_open_read(wfk0,wfkfile,formeig0,iomode0,get_unit(),xmpi_comm_self) 237 ! Get header from the gs file 238 call hdr_copy(wfk0%hdr, hdr) 239 240 ! Identify the type of RF Wavefunction files 241 use_ncevk = .False. 242 do ii=1,3 243 use_ncevk(ii) = endswith(infiles(ii), "_EVK.nc") 244 end do 245 246 ! Read ddk here from WFK files or from EVK.nc (only the header in the latter case) 247 do ii=1,3 248 249 call nctk_fort_or_ncfile(infiles(ii), iomode_ddk(ii), msg) 250 if (len_trim(msg) /= 0) ABI_ERROR(msg) 251 if (iomode_ddk(ii) == IO_MODE_MPI) iomode_ddk(ii) = IO_MODE_FORTRAN 252 253 if (.not. use_ncevk(ii)) then 254 call wfk_open_read(wfks(ii), infiles(ii), formeig1, iomode_ddk(ii), get_unit(), xmpi_comm_self) 255 call hdr_copy(wfks(ii)%hdr, hdr_ddk(ii)) 256 else 257 258 #ifdef HAVE_NETCDF 259 NCF_CHECK(nctk_open_read(ncid, infiles(ii), xmpi_comm_self)) 260 call hdr_ncread(hdr_ddk(ii), ncid, fform) 261 ABI_CHECK(fform /= 0, sjoin("Error while reading:", infiles(ii))) 262 263 NCF_CHECK(nf90_close(ncid)) 264 #else 265 ABI_ERROR("Netcdf not available!") 266 #endif 267 268 end if 269 end do 270 271 ! if(any(iomode_ddk(:)/=iomode0))then 272 ! write(msg, "(5a)")& 273 !& ' The ground-state and ddk files should have the same format,',ch10,& 274 !& ' either FORTRAN binary or NetCDF, which is not the case.',ch10,& 275 !& ' Action : see input variable iomode.' 276 ! ABI_ERROR(msg) 277 ! endif 278 279 ! Perform basic consistency tests for the GS WFK and the DDK files, e.g. 280 ! k-points and their order, spins, number of bands could differ in the four files. 281 ! Note indeed that we must have the same quantities in all the files. 282 283 if (.not. use_ncevk(1)) then 284 285 write(msg, "(12a)")ch10,& 286 ' Check the consistency of the wavefunction files (esp. k point and number of bands). ',ch10,& 287 ' Will compare, pairwise ( 1/2, 2/3, 3/4 ), the four following files :',ch10,trim(wfkfile) 288 ! split the write since long filenames can bust the 500 char limit of 'msg' 289 call wrtout(std_out, msg) 290 do ii=1,3 291 write(msg, "(12a)")trim(infiles(ii)) 292 call wrtout(std_out, msg) 293 enddo 294 295 if (hdr%compare(hdr_ddk(1)) /= 0) then 296 write(msg, "(3a)")" GS WFK file and ddkfile ",trim(infiles(1))," are not consistent. See above messages." 297 ABI_ERROR(msg) 298 end if 299 do ii=1,2 300 if (wfks(ii)%compare(wfks(ii+1)) /= 0) then 301 write(msg, "(2(a,i0,a))")" ddkfile", ii," and ddkfile ",ii+1, ", are not consistent. See above messages" 302 ABI_ERROR(msg) 303 end if 304 enddo 305 endif 306 307 ! TODO: one should perform basic consistency tests for the EVK files, e.g. 308 ! k-points and their order, spins, number of bands could differ in the four files. 309 ! Note indeed that we are assuming the same numer of bands in all the files. 310 311 !Handle electron-phonon file 312 ep_nc_fname = 'test_EP.nc'; if (do_temperature) ep_nc_fname = epfile 313 do_ep_renorm = file_exists(ep_nc_fname) 314 ep_ntemp = 1 315 if (do_ep_renorm) then 316 call eprenorms_from_epnc(Epren,ep_nc_fname) 317 ep_ntemp = Epren%ntemp 318 else if (do_temperature) then 319 ABI_ERROR("You have asked for temperature but the epfile is not present !") 320 end if 321 322 ! autoparal section 323 if (autoparal /= 0 .and. max_ncpus /= 0) then 324 write(std_out,'(a)')"--- !Autoparal" 325 write(std_out,"(a)")'#Autoparal section for Optic runs.' 326 write(std_out,"(a)") "info:" 327 write(std_out,"(a,i0)")" autoparal: ",autoparal 328 write(std_out,"(a,i0)")" max_ncpus: ",max_ncpus 329 write(std_out,"(a,i0)")" nspinor: ",hdr%nspinor 330 write(std_out,"(a,i0)")" nsppol: ",hdr%nsppol 331 write(std_out,"(a,i0)")" nkpt: ",hdr%nkpt 332 write(std_out,"(a,i0)")" mband: ",maxval(hdr%nband) 333 334 work_size = hdr%nkpt !* hdr%nsppol 335 336 ! List of configurations. 337 write(std_out,"(a)")"configurations:" 338 do ii=1,max_ncpus 339 if (ii > work_size) cycle 340 nks_per_proc = work_size / ii 341 nks_per_proc = nks_per_proc + MOD(work_size, ii) 342 eff = (one * work_size) / (ii * nks_per_proc) 343 write(std_out,"(a,i0)")" - tot_ncpus: ",ii !* omp_ncpus 344 write(std_out,"(a,i0)")" mpi_ncpus: ",ii 345 write(std_out,"(a,i0)")" omp_ncpus: ",1 346 write(std_out,"(a,f12.9)")" efficiency: ",eff 347 !write(,"(a,f12.2)")" mem_per_cpu: ",mempercpu_mb 348 end do 349 350 write(std_out,'(a)')"..." 351 ABI_ERROR_NODUMP("aborting now") 352 end if 353 354 end if ! my_rank == master 355 356 ! Master broadcasts input variables. 357 call hdr%bcast(master, my_rank, comm) 358 call xmpi_bcast(broadening, master, comm, ierr) 359 call xmpi_bcast(domega, master, comm, ierr) 360 call xmpi_bcast(maxomega, master, comm, ierr) 361 call xmpi_bcast(scissor, master, comm, ierr) 362 call xmpi_bcast(tolerance, master, comm, ierr) 363 call xmpi_bcast(num_lin_comp, master, comm, ierr) 364 call xmpi_bcast(prtlincompmatrixelements, master, comm, ierr) 365 call xmpi_bcast(nband_sum, master, comm, ierr) 366 call xmpi_bcast(lin_comp,master, comm, ierr) 367 call xmpi_bcast(num_nonlin_comp,master, comm, ierr) 368 call xmpi_bcast(nonlin_comp, master, comm, ierr) 369 call xmpi_bcast(num_linel_comp, master, comm, ierr) 370 call xmpi_bcast(linel_comp, master, comm, ierr) 371 call xmpi_bcast(num_nonlin2_comp, master, comm, ierr) 372 call xmpi_bcast(nonlin2_comp, master, comm, ierr) 373 call xmpi_bcast(do_antiresonant, master, comm, ierr) 374 call xmpi_bcast(do_ep_renorm, master, comm, ierr) 375 call xmpi_bcast(ep_ntemp, master, comm, ierr) 376 call xmpi_bcast(filnam_out, master, comm, ierr) 377 call xmpi_bcast(prefix, master, comm, ierr) 378 if (do_ep_renorm) call eprenorms_bcast(Epren, master, comm) 379 380 ! Extract basic info from the header 381 bantot = hdr%bantot 382 nkpt = hdr%nkpt 383 nsppol = hdr%nsppol 384 385 ! Get mband as the maximum value of nband(nkpt) and init nband_sum if negative 386 mband = maxval(hdr%nband) 387 ABI_CHECK(all(hdr%nband == mband), "nband must be constant across kpts") 388 if (nband_sum == -1) nband_sum = mband 389 if (nband_sum <= 0 .or. nband_sum > mband) then 390 ABI_ERROR(sjoin("nband_sum should be in [1, mband] while it is:", itoa(nband_sum), "with mband:", itoa(mband))) 391 end if 392 393 ! Initializes crystal object 394 call crystal_init(hdr%amu, cryst, 0, hdr%natom, hdr%npsp, hdr%ntypat, & 395 hdr%nsym, hdr%rprimd, hdr%typat, hdr%xred, hdr%zionpsp, hdr%znuclpsp, 1, & 396 (hdr%nspden==2 .and. hdr%nsppol==1),remove_inv, hdr%title,& 397 hdr%symrel, hdr%tnons, hdr%symafm) 398 399 if (my_rank == master) then 400 write(std_out,*) 401 write(std_out,'(a,3f10.5,a)' )' rprimd(bohr) =',cryst%rprimd(1:3,1) 402 write(std_out,'(a,3f10.5,a)' )' ',cryst%rprimd(1:3,2) 403 write(std_out,'(a,3f10.5,a)' )' ',cryst%rprimd(1:3,3) 404 write(std_out,'(a,i8)') ' natom =',cryst%natom 405 write(std_out,'(a,2i8)') ' nkpt,mband =',nkpt,mband 406 write(std_out,'(a, f10.5,a)' ) ' ecut =',hdr%ecut_eff,' Ha' 407 end if 408 409 ! Read the eigenvalues of ground-state and ddk files 410 ABI_MALLOC(eigen0, (mband*nkpt*nsppol)) 411 ! MG: Do not understand why not [...,3] 412 ABI_MALLOC(eigen11, (2*mband*mband*nkpt*nsppol)) 413 ABI_MALLOC(eigen12, (2*mband*mband*nkpt*nsppol)) 414 ABI_MALLOC(eigen13, (2*mband*mband*nkpt*nsppol)) 415 416 if (my_rank == master) then 417 ABI_MALLOC(eigtmp, (2*mband*mband)) 418 ABI_MALLOC(eig0tmp, (mband)) 419 420 do ii=1,3 421 if (.not. use_ncevk(ii)) cycle 422 #ifdef HAVE_NETCDF 423 NCF_CHECK(nctk_open_read(ncid, infiles(ii), xmpi_comm_self)) 424 varid = nctk_idname(ncid, "h1_matrix_elements") 425 outeig => eigen11 426 if (ii == 2) outeig => eigen12 427 if (ii == 3) outeig => eigen13 428 NCF_CHECK(nf90_get_var(ncid, varid, outeig, count=[2, mband, mband, nkpt, nsppol])) 429 NCF_CHECK(nf90_close(ncid)) 430 #else 431 ABI_ERROR("Netcdf not available!") 432 #endif 433 end do 434 435 bdtot0_index=0 ; bdtot_index=0 436 do isppol=1,nsppol 437 do ikpt=1,nkpt 438 nband1 = hdr%nband(ikpt+(isppol-1)*nkpt) 439 eigtmp = zero 440 eig0tmp = zero 441 442 call wfk0%read_eigk(ikpt,isppol,xmpio_single,eig0tmp) 443 eigen0(1+bdtot0_index:nband1+bdtot0_index)=eig0tmp(1:nband1) 444 445 ! Read DDK matrix elements from WFK 446 do ii=1,3 447 if (.not. use_ncevk(ii)) then 448 call wfks(ii)%read_eigk(ikpt, isppol, xmpio_single, eigtmp) 449 if (ii == 1) eigen11(1+bdtot_index:2*nband1**2+bdtot_index)=eigtmp(1:2*nband1**2) 450 if (ii == 2) eigen12(1+bdtot_index:2*nband1**2+bdtot_index)=eigtmp(1:2*nband1**2) 451 if (ii == 3) eigen13(1+bdtot_index:2*nband1**2+bdtot_index)=eigtmp(1:2*nband1**2) 452 !ABI_CHECK(wfks(ii)%nband(ikpt,isppol) == nband1, "ddk1 nband1") 453 end if 454 end do 455 bdtot0_index=bdtot0_index+nband1 456 bdtot_index=bdtot_index+2*nband1**2 457 end do 458 end do 459 460 call wfk0%close() 461 do ii=1,3 462 if (.not. use_ncevk(ii)) call wfks(ii)%close() 463 end do 464 465 ABI_FREE(eigtmp) 466 ABI_FREE(eig0tmp) 467 end if ! master 468 469 call xmpi_bcast(eigen0, master,comm, ierr) 470 call xmpi_bcast(eigen11, master, comm, ierr) 471 call xmpi_bcast(eigen12, master, comm, ierr) 472 call xmpi_bcast(eigen13, master, comm, ierr) 473 474 ! Recompute fermie from header 475 ! WARNING no guarantee that it works for other materials than insulators 476 477 ABI_MALLOC(doccde, (mband * nkpt * nsppol)) 478 479 call ebands_init(bantot, ks_ebands, hdr%nelect, hdr%ne_qFD, hdr%nh_qFD, hdr%ivalence,& 480 doccde, eigen0, hdr%istwfk, hdr%kptns, & 481 hdr%nband, nkpt, hdr%npwarr, nsppol, hdr%nspinor, hdr%tphysel, broadening, hdr%occopt, hdr%occ, hdr%wtk, & 482 hdr%cellcharge, hdr%kptopt, hdr%kptrlatt_orig, hdr%nshiftk_orig, hdr%shiftk_orig, & 483 hdr%kptrlatt, hdr%nshiftk, hdr%shiftk) 484 485 ABI_FREE(eigen0) 486 ABI_FREE(doccde) 487 !ks_ebands = ebands_from_hdr(hdr, mband, ene3d, nelect) result(ebands) 488 489 !YG : should we use broadening for ebands_init 490 call ebands_update_occ(ks_ebands, -99.99d0) 491 492 !size of the frequency range 493 nomega=int((maxomega+domega*0.001_dp)/domega) 494 maxomega = dble(nomega)*domega 495 496 optic_ncid = nctk_noid 497 if (my_rank == master) then 498 write(std_out,'(a,f10.5,a,f10.5,a)' )' fermie =',ks_ebands%fermie,' Ha',ks_ebands%fermie*Ha_eV,' eV' 499 write(std_out,'(a,f10.5,a)')' Scissor shift =', scissor, ' Ha' 500 write(std_out,'(a,f10.5,a)')' Tolerance on closeness to singularities =', tolerance, ' Ha' 501 write(std_out,'(a,f10.5,a)')' Smearing factor =', broadening, ' Ha' 502 if (do_antiresonant) then 503 write(std_out,'(a)') ' Will use the antiresonant approximation ' 504 else 505 write(std_out,'(a)') ' Will not use the antiresonant approximation (only available for nonlin2 and linel components!) ' 506 end if 507 write(std_out,'(a)') ' linear coeffs to be calculated : ' 508 write(std_out,'(9i3)') lin_comp(1:num_lin_comp) 509 write(std_out,'(a)') ' non-linear coeffs to be calculated : ' 510 write(std_out,'(27i4)') nonlin_comp(1:num_nonlin_comp) 511 write(std_out,'(a)') ' electronic part of electro-optic coeffs to be calculated :' 512 write(std_out,'(27i4)') linel_comp(1:num_linel_comp) 513 write(std_out,'(a)') ' non-linear coeffs (V2) to be calculated :' 514 write(std_out,'(27i4)') nonlin2_comp(1:num_nonlin2_comp) 515 write(std_out,'(a,i1)') ' linear optic matrix elements will be printed :',prtlincompmatrixelements 516 517 #ifdef HAVE_NETCDF 518 ! Open netcdf file that will contain output results (only master is supposed to write) 519 NCF_CHECK_MSG(nctk_open_create(optic_ncid, strcat(prefix, "_OPTIC.nc"), xmpi_comm_self), "Creating _OPTIC.nc") 520 521 ! Add header, crystal, and ks_ebands 522 ! Note that we write the KS bands without EPH interaction (if any). 523 NCF_CHECK(hdr%ncwrite(optic_ncid, 666, nc_define=.True.)) 524 NCF_CHECK(cryst%ncwrite(optic_ncid)) 525 NCF_CHECK(ebands_ncwrite(ks_ebands, optic_ncid)) 526 527 ! Add optic input variables. 528 NCF_CHECK(nctk_def_dims(optic_ncid, [nctkdim_t("ntemp", ep_ntemp), nctkdim_t("nomega", nomega)], defmode=.True.)) 529 530 ncerr = nctk_def_iscalars(optic_ncid, [character(len=nctk_slen) :: & 531 "do_antiresonant", "do_ep_renorm", "nband_sum"]) 532 NCF_CHECK(ncerr) 533 ncerr = nctk_def_dpscalars(optic_ncid, [character(len=nctk_slen) :: & 534 "broadening", "domega", "maxomega", "scissor", "tolerance"]) 535 NCF_CHECK(ncerr) 536 537 ! Define arrays containing output results 538 ncerr = nctk_def_arrays(optic_ncid, [nctkarr_t('wmesh', "dp", "nomega")]) 539 NCF_CHECK(ncerr) 540 541 if (num_lin_comp > 0) then 542 ! Linear optic results. 543 NCF_CHECK(nctk_def_dims(optic_ncid, nctkdim_t("linopt_ncomp", num_lin_comp))) 544 ncerr = nctk_def_arrays(optic_ncid, [ & 545 nctkarr_t('linopt_components', "int", "linopt_ncomp"), & 546 nctkarr_t('linopt_epsilon', "dp", "two, nomega, linopt_ncomp, ntemp") & 547 ]) 548 NCF_CHECK(ncerr) 549 if (prtlincompmatrixelements == 1) then 550 ! Linear optic matrix elements 551 ncerr = nctk_def_dims(optic_ncid, [ & 552 nctkdim_t("nkpt", nkpt), & 553 nctkdim_t("nband", mband), & 554 nctkdim_t("nsppol", nsppol)], & 555 defmode=.True.) 556 NCF_CHECK(ncerr) 557 ncerr = nctk_def_arrays(optic_ncid, [ & 558 !nctkarr_t('linopt_components', "int", "linopt_ncomp"), & 559 nctkarr_t('linopt_matrix_elements', "dp", "two, nband, nband, nkpt, nsppol, linopt_ncomp, ntemp"), & 560 nctkarr_t('linopt_renorm_eigs', "dp", "two, nband, nkpt, nsppol"), & 561 nctkarr_t('linopt_occupations', "dp", "nband, nkpt, nsppol"), & 562 nctkarr_t('linopt_wkpts', "dp", "nkpt") & 563 ]) 564 NCF_CHECK(ncerr) 565 endif 566 end if 567 568 if (num_nonlin_comp > 0) then 569 ! Second harmonic generation. 570 NCF_CHECK(nctk_def_dims(optic_ncid, nctkdim_t("shg_ncomp", num_nonlin_comp))) 571 ncerr = nctk_def_arrays(optic_ncid, [ & 572 nctkarr_t('shg_components', "int", "shg_ncomp"), & 573 nctkarr_t('shg_inter2w', "dp", "two, nomega, shg_ncomp, ntemp"), & 574 nctkarr_t('shg_inter1w', "dp", "two, nomega, shg_ncomp, ntemp"), & 575 nctkarr_t('shg_intra2w', "dp", "two, nomega, shg_ncomp, ntemp"), & 576 nctkarr_t('shg_intra1w', "dp", "two, nomega, shg_ncomp, ntemp"), & 577 nctkarr_t('shg_intra1wS', "dp", "two, nomega, shg_ncomp, ntemp"), & 578 nctkarr_t('shg_chi2tot', "dp", "two, nomega, shg_ncomp, ntemp") & 579 ]) 580 NCF_CHECK(ncerr) 581 end if 582 583 if (num_linel_comp > 0) then 584 ! linear electro-optic (LEO) susceptibility 585 NCF_CHECK(nctk_def_dims(optic_ncid, nctkdim_t("leo_ncomp", num_linel_comp))) 586 ncerr = nctk_def_arrays(optic_ncid, [ & 587 nctkarr_t('leo_components', "int", "leo_ncomp"), & 588 nctkarr_t('leo_chi', "dp", "two, nomega, leo_ncomp, ntemp"), & 589 nctkarr_t('leo_eta', "dp", "two, nomega, leo_ncomp, ntemp"), & 590 nctkarr_t('leo_sigma', "dp", "two, nomega, leo_ncomp, ntemp"), & 591 nctkarr_t('leo_chi2tot', "dp", "two, nomega, leo_ncomp, ntemp") & 592 ]) 593 NCF_CHECK(ncerr) 594 end if 595 596 if (num_nonlin2_comp > 0) then 597 ! non-linear electro-optic susceptibility 598 NCF_CHECK(nctk_def_dims(optic_ncid, nctkdim_t("leo2_ncomp", num_nonlin2_comp))) 599 ncerr = nctk_def_arrays(optic_ncid, [ & 600 nctkarr_t('leo2_components', "int", "leo2_ncomp"), & 601 nctkarr_t('leo2_chiw', "dp", "two, nomega, leo2_ncomp, ntemp"), & 602 nctkarr_t('leo2_etaw', "dp", "two, nomega, leo2_ncomp, ntemp"), & 603 nctkarr_t('leo2_chi2w', "dp", "two, nomega, leo2_ncomp, ntemp"), & 604 nctkarr_t('leo2_eta2w', "dp", "two, nomega, leo2_ncomp, ntemp"), & 605 nctkarr_t('leo2_sigmaw', "dp", "two, nomega, leo2_ncomp, ntemp"), & 606 nctkarr_t('leo2_chi2tot', "dp", "two, nomega, leo2_ncomp, ntemp") & 607 ]) 608 NCF_CHECK(ncerr) 609 end if 610 611 NCF_CHECK(nctk_set_datamode(optic_ncid)) 612 613 ! Write wmesh here. 614 ABI_MALLOC(wmesh, (nomega)) 615 do ii=1,nomega 616 ! This to be consistent with the value used in m_optic_tools 617 ! In principle wmesh should be passed to the children and a lot of code 618 ! should be rewritten to be more cache-friendly ... 619 wmesh(ii) = (ii-1)*domega * (13.60569172*2._dp) 620 end do 621 NCF_CHECK(nf90_put_var(optic_ncid, nctk_idname(optic_ncid, "wmesh"), wmesh)) 622 ABI_FREE(wmesh) 623 624 if (num_lin_comp > 0) then 625 NCF_CHECK(nf90_put_var(optic_ncid, nctk_idname(optic_ncid, "linopt_components"), lin_comp(1:num_lin_comp))) 626 end if 627 if (num_nonlin_comp > 0) then 628 NCF_CHECK(nf90_put_var(optic_ncid, nctk_idname(optic_ncid, "shg_components"), nonlin_comp(1:num_nonlin_comp))) 629 end if 630 if (num_linel_comp > 0) then 631 NCF_CHECK(nf90_put_var(optic_ncid, nctk_idname(optic_ncid, "leo_components"), linel_comp(1:num_linel_comp))) 632 end if 633 if (num_nonlin2_comp > 0) then 634 NCF_CHECK(nf90_put_var(optic_ncid, nctk_idname(optic_ncid, "leo2_components"), nonlin2_comp(1:num_nonlin2_comp))) 635 end if 636 637 ! Write optic input variables. 638 ii = 0; if (do_antiresonant) ii = 1 639 jj = 0; if (do_ep_renorm) jj = 1 640 ncerr = nctk_write_iscalars(optic_ncid, [character(len=nctk_slen) :: & 641 "do_antiresonant", "do_ep_renorm", "nband_sum"], & 642 [ii, jj, nband_sum]) 643 NCF_CHECK(ncerr) 644 645 ncerr = nctk_write_dpscalars(optic_ncid, [character(len=nctk_slen) :: & 646 "broadening", "domega", "maxomega", "scissor", "tolerance"], & 647 [broadening, domega, maxomega, scissor, tolerance]) 648 NCF_CHECK(ncerr) 649 #endif 650 end if 651 652 ! Get velocity matrix elements in cartesian coordinates from reduced coords. 653 call wrtout(std_out," optic : Call pmat2cart") 654 ABI_MALLOC(pmat, (mband, mband, nkpt, 3, nsppol)) 655 call pmat2cart(eigen11, eigen12, eigen13, mband, nkpt, nsppol, pmat, cryst%rprimd) 656 ABI_FREE(eigen11) 657 ABI_FREE(eigen12) 658 ABI_FREE(eigen13) 659 660 ! Renormalize matrix elements if scissors is being used. 661 call pmat_renorm(ks_ebands%fermie, ks_ebands%eig, mband, nkpt, nsppol, pmat, scissor) 662 663 !--------------------------------------------------------------------------------- 664 ! Perform calculations 665 !--------------------------------------------------------------------------------- 666 667 ! XG_2020_05_25 : All these subroutines should be rationalized. There are numerous 668 ! similar sections, e.g. at the level of the checking, and set up ... 669 670 ! v1,v2=desired component of the dielectric function(integer) 1=x,2=y,3=z 671 ! nmesh=desired number of energy mesh points(integer) 672 ! de=desired step in energy(real); nmesh*de=maximum energy 673 ! scissor=scissors shift in Ha(real) 674 ! brod=broadening in Ha(real) 675 676 ! optical frequency dependent dielectric function for semiconductors 677 call wrtout(std_out," optic : Call linopt") 678 679 do itemp=1,ep_ntemp 680 call ebands_copy(ks_ebands, eph_ebands) 681 if (do_ep_renorm) call renorm_bst(Epren, eph_ebands, cryst, itemp, do_lifetime=.True.,do_check=.True.) 682 do ii=1,num_lin_comp 683 lin1 = int(lin_comp(ii)/10.0_dp) 684 lin2 = mod(lin_comp(ii),10) 685 write(msg,*) ' linopt ', lin1,lin2 686 call wrtout(std_out, msg) 687 call int2char4(lin1,s1) 688 call int2char4(lin2,s2) 689 call int2char4(itemp,stemp) 690 ABI_CHECK((s1(1:1)/='#'),'Bug: string length too short!') 691 ABI_CHECK((s2(1:1)/='#'),'Bug: string length too short!') 692 ABI_CHECK((stemp(1:1)/='#'),'Bug: string length too short!') 693 tmp_radix = trim(prefix)//"_"//trim(s1)//"_"//trim(s2) 694 if (do_ep_renorm) tmp_radix = trim(prefix)//"_"//trim(s1)//"_"//trim(s2)//"_T"//trim(stemp) 695 call linopt(ii, itemp, nband_sum, cryst, ks_ebands, eph_ebands, pmat, & 696 lin1, lin2, nomega, domega, scissor, broadening, tmp_radix, optic_ncid, prtlincompmatrixelements, comm) 697 end do 698 call ebands_free(eph_ebands) 699 end do 700 701 if (do_ep_renorm) call eprenorms_free(Epren) 702 703 ! second harmonic generation susceptibility for semiconductors 704 call wrtout(std_out," optic : Call nlinopt") 705 do ii=1,num_nonlin_comp 706 nlin1 = int( nonlin_comp(ii)/100.0_dp) 707 nlin2 = int((nonlin_comp(ii)-nlin1*100.0_dp)/10.0_dp) 708 nlin3 = mod( nonlin_comp(ii),10) 709 write(msg,*) ' nlinopt ', nlin1,nlin2,nlin3 710 call wrtout(std_out, msg) 711 call int2char4(nlin1,s1) 712 call int2char4(nlin2,s2) 713 call int2char4(nlin3,s3) 714 ABI_CHECK((s1(1:1)/='#'),'Bug: string length too short!') 715 ABI_CHECK((s2(1:1)/='#'),'Bug: string length too short!') 716 ABI_CHECK((s3(1:1)/='#'),'Bug: string length too short!') 717 tmp_radix = trim(prefix)//"_"//trim(s1)//"_"//trim(s2)//"_"//trim(s3) 718 itemp = 1 719 720 if (hdr%kptopt == 1) then 721 ABI_WARNING("second harmonic generation with symmetries (kptopt == 1) is not tested. Use at your own risk!") 722 end if 723 724 call nlinopt(ii, itemp, nband_sum, cryst, ks_ebands, pmat, & 725 nlin1, nlin2, nlin3, nomega, domega, scissor, broadening, tolerance, tmp_radix, optic_ncid, comm) 726 end do 727 728 ! linear electro-optic susceptibility for semiconductors 729 call wrtout(std_out," optic : Call linelop") 730 do ii=1,num_linel_comp 731 linel1 = int(linel_comp(ii)/100.0_dp) 732 linel2 = int((linel_comp(ii)-linel1*100.0_dp)/10.0_dp) 733 linel3 = mod(linel_comp(ii),10) 734 write(msg,*) ' linelop ',linel1,linel2,linel3 735 call wrtout(std_out, msg) 736 call int2char4(linel1,s1) 737 call int2char4(linel2,s2) 738 call int2char4(linel3,s3) 739 tmp_radix = trim(prefix)//"_"//trim(s1)//"_"//trim(s2)//"_"//trim(s3) 740 itemp = 1 741 742 if (hdr%kptopt == 1) then 743 ABI_ERROR("linear electro-optic with symmetries (kptopt == 1) is not tested. Use at your own risk!") 744 end if 745 746 call linelop(ii, itemp, nband_sum, cryst, ks_ebands, pmat, & 747 linel1, linel2, linel3, nomega, domega, scissor, broadening, & 748 tolerance, tmp_radix, do_antiresonant, optic_ncid, comm) 749 end do 750 751 ! nonlinear electro-optical susceptibility for semiconductors 752 call wrtout(std_out," optic : Call nonlinopt") 753 do ii=1,num_nonlin2_comp 754 nonlin1 = int(nonlin2_comp(ii)/100.0_dp) 755 nonlin2 = int((nonlin2_comp(ii)-nonlin1*100.0_dp)/10.0_dp) 756 nonlin3 = mod(nonlin2_comp(ii),10) 757 write(msg,*) ' nonlinopt ',nonlin1,nonlin2,nonlin3 758 call wrtout(std_out, msg) 759 call int2char4(nonlin1,s1) 760 call int2char4(nonlin2,s2) 761 call int2char4(nonlin3,s3) 762 tmp_radix = trim(prefix)//"_"//trim(s1)//"_"//trim(s2)//"_"//trim(s3) 763 itemp = 1 764 765 if (hdr%kptopt == 1) then 766 ABI_ERROR("nonlinear electro-optic with symmetries (kptopt == 1) is not tested. Use at your own risk!") 767 end if 768 769 call nonlinopt(ii, itemp, nband_sum, cryst, ks_ebands, pmat, & 770 nonlin1, nonlin2, nonlin3, nomega, domega, scissor, broadening, tolerance, tmp_radix, & 771 do_antiresonant, optic_ncid, comm) 772 end do 773 774 ! Free memory 775 ABI_FREE(pmat) 776 call hdr%free() 777 do ii=1,3 778 call hdr_ddk(ii)%free() 779 end do 780 call ebands_free(ks_ebands) 781 call cryst%free() 782 783 call timein(tcpu, twall) 784 785 tsec(1) = tcpu - tcpui 786 tsec(2) = twall - twalli 787 788 if (my_rank == master) then 789 write(std_out,'(a,80a,a,a,a)' )ch10,('=',ii=1,80),ch10,ch10,' Calculation completed.' 790 write(std_out, '(a,a,a,f13.1,a,f13.1)' ) & 791 '-',ch10,'- Proc. 0 individual time (sec): cpu=',tsec(1),' wall=',tsec(2) 792 end if 793 794 call xmpi_sum(tsec, comm, ierr) 795 796 if (my_rank == master) then 797 ! Write YAML document with the final summary. 798 ! we use this doc to test whether the calculation is completed. 799 write(std_out,"(a)")"--- !FinalSummary" 800 write(std_out,"(a)")"program: optic" 801 write(std_out,"(2a)")"version: ",trim(abinit_version) 802 write(std_out,"(2a)")"start_datetime: ",start_datetime 803 write(std_out,"(2a)")"end_datetime: ",asctime() 804 write(std_out,"(a,f13.1)")"overall_cpu_time: ",tsec(1) 805 write(std_out,"(a,f13.1)")"overall_wall_time: ",tsec(2) 806 write(std_out,"(a,i0)")"mpi_procs: ",xmpi_comm_size(xmpi_world) 807 write(std_out,"(a,i0)")"omp_threads: ",xomp_get_num_threads(open_parallel=.True.) 808 !write(std_out,"(a,i0)")"num_warnings: ",nwarning 809 !write(std_out,"(a,i0)")"num_comments: ",ncomment 810 write(std_out,"(a)")"..." 811 call flush_unit(std_out) 812 end if 813 814 #ifdef HAVE_NETCDF 815 if (my_rank == master) then 816 NCF_CHECK(nf90_close(optic_ncid)) 817 end if 818 #endif 819 820 ! Write information on file about the memory before ending mpi module, if memory profiling is enabled 821 call abinit_doctor(filnam) 822 823 100 call xmpi_end() 824 825 end program optic