TABLE OF CONTENTS


ABINIT/optic [ Programs ]

[ Top ] [ 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-2022 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_build_info
 49  use m_optic_tools
 50  use m_wfk
 51  use m_nctk
 52  use m_hdr
 53  use m_ebands
 54  use m_eprenorms
 55  use m_crystal
 56  use m_argparse
 57 #ifdef HAVE_NETCDF
 58  use netcdf
 59 #endif
 60 
 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