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 (C) 2002-2018 ABINIT group (SSharma,MVer,VRecoules,YG)
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

```  bantot
doccde(mband*nkpt_rbz*nsppol)=derivative of occ_rbz wrt the energy.
domega=frequency range
eigen0(mband*nkpt_rbz*nsppol)=GS eigenvalues at k (hartree).
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
ecut=kinetic energy planewave cutoff (hartree).
entropy= entropy associated with the smearing (adimensional)
fermie= fermi energy (Hartree)
gmet(3,3)=reciprocal space metric (\$\textrm{bohr}^{2}\$).
gmet_inv(3,3)=inverse of reciprocal space metric (\$\textrm{bohr}^{2}\$).
gprimd(3,3)=dimensional primitive translations for reciprocal space(bohr^-1).
nomega=number of frequency for conductivity computation
mband=maximum number of bands.
natom = number of atoms in the unit cell.
nband(nkpt*nsppol)=number of bands at each RF k point for each spin.
nelect=number of electrons per unit cell
nkpt=number of k points in the IBZ for this perturbation
ngfft(3)=integer fft box dimensions.
nspinor=number of spinorial components of the wavefunctions.
nsppol=1 for unpolarized, 2 for spin-polarized.
ntypat = number of atom types.
occ(mband*nkpt*nsppol)=occupation number for each band and k.
occopt==option for occupancies
rmet(3,3)=real space metric (\$\textrm{bohr}^{2}\$).
rprimd(3,3)=real space primitive translations.
of primitive translations.
broadening=smearing width (or temperature) in Hartree
ucvol=unit cell volume in (\$\textrm{bohr}^{3}\$).
maxomega=frequency windows for computations of sigma
wtk(nkpt)=weight assigned to each k point.
znucl(natom)=atomic number of atoms
```

PARENTS

CHILDREN

```      abi_io_redirect,abimem_init,abinit_doctor,crystal_free,crystal_init
ebands_copy,ebands_free,ebands_init,ebands_update_occ,eprenorms_bcast
eprenorms_free,eprenorms_from_epnc,flush_unit,hdr_bcast,hdr_copy
hdr_free,herald,int2char4,linelop,linopt,mati3inv,matr3inv,metric
nctk_fort_or_ncfile,nlinopt,nonlinopt,pmat2cart,pmat_renorm,renorm_bst
xmpi_end,xmpi_init,xmpi_sum
```

SOURCE

``` 74 #if defined HAVE_CONFIG_H
75 #include "config.h"
76 #endif
77
78 #include "abi_common.h"
79
80 program optic
81
82  use defs_basis
83  use defs_datatypes
84  use defs_abitypes
85  use m_errors
86  use m_xmpi
87  use m_xomp
88  use m_abicore
89  use m_build_info
90  use m_optic_tools
91  use m_wfk
92  use m_nctk
93  use m_hdr
94  use m_ebands
95  use m_eprenorms
96  use m_crystal
97 #ifdef HAVE_NETCDF
98  use netcdf
99 #endif
100
101  use m_specialmsg,     only : specialmsg_getcount, herald
102  use m_time ,          only : asctime, timein
103  use m_symtk,          only : mati3inv, matr3inv
104  use m_geometry,       only : metric
105  use m_io_tools,       only : flush_unit, open_file, file_exists, get_unit
106  use m_numeric_tools,  only : c2r
107  use m_fstrings,       only : int2char4, itoa, sjoin, strcat, endswith
108  use m_crystal_io,     only : crystal_ncwrite
109
110 !This section has been created automatically by the script Abilint (TD).
111 !Do not modify the following lines by hand.
112 #undef ABI_FUNC
113 #define ABI_FUNC 'optic'
114 !End of the abilint section
115
116  implicit none
117
118 !Arguments -----------------------------------
119
120 !Local variables-------------------------------
121  integer,parameter :: formeig0=0,formeig1=1,tim_rwwf=0,master=0
122  integer :: finunt,ep_ntemp,itemp
123  integer :: bantot,bdtot0_index,bdtot_index
125  integer :: isppol,mband,nomega,natom,nband1,nsym
126  integer :: nkpt,nspinor,nsppol,ntypat
127  integer :: occopt,nks_per_proc,work_size,lin1,lin2,nlin1,nlin2,nlin3
128  integer :: linel1,linel2,linel3,nonlin1,nonlin2,nonlin3
129  integer :: iomode,comm,nproc,my_rank, optic_ncid
130 #ifdef HAVE_NETCDF
131  integer :: ncid, varid, ncerr
132 #endif
133  integer :: num_lin_comp=1,num_nonlin_comp=0,num_linel_comp=0,num_nonlin2_comp=0
134  integer :: autoparal=0,max_ncpus=0
135  integer :: nonlin_comp(27) = 0, linel_comp(27) = 0, nonlin2_comp(27) = 0
136  integer :: lin_comp(9) = [11, 22 ,33, 12, 13, 21, 23, 31, 32]
137  real(dp) :: domega,ecut,fermie, eff
139  real(dp) :: tcpu,tcpui,twall,twalli,nelect
140  logical :: do_antiresonant, do_temperature
141  logical :: do_ep_renorm
142  logical,parameter :: remove_inv = .False.
143  type(hdr_type) :: hdr
144  type(ebands_t) :: ks_ebands, eph_ebands
145  type(crystal_t) :: cryst
146  type(eprenorms_t) :: Epren
147 !arrays
148  integer,allocatable :: istwfk(:), npwarr(:), nband(:), symrel(:,:,:), symrec(:,:,:)
149  real(dp) :: tsec(2), gmet(3,3),gmet_inv(3,3),gprimd(3,3),rmet(3,3),rprimd(3,3),gprimd_trans(3,3)
150  real(dp),allocatable :: symcart(:,:,:), kpt(:,:),wmesh(:)
151  real(dp),allocatable :: cond_kg(:),cond_nd(:),doccde(:), eig0tmp(:), eigen0(:)
152  real(dp),target,allocatable :: eigen11(:),eigen12(:),eigen13(:)
153  real(dp),allocatable :: occ(:),wtk(:),eigtmp(:)
154  real(dp), ABI_CONTIGUOUS pointer :: outeig(:)
155  complex(dpc),allocatable :: pmat(:,:,:,:,:)
156  logical :: use_ncddk(0:3)
157  character(len=fnlen) :: filnam,wfkfile,ddkfile_1,ddkfile_2,ddkfile_3,filnam_out, epfile
158  character(len=fnlen) :: infiles(0:3)
159 ! for the moment this is imposed by the format in linopt.f and nlinopt.f
161  character(len=10) :: s1,s2,s3,stemp
162  character(len=24) :: codename
163  character(len=24) :: start_datetime
164  character(len=500) :: msg
165  character(len=fnlen) :: ep_nc_fname
166  type(wfk_t) :: wfk0
167  type(wfk_t) :: wfks(0:3)
168
169  ! Input file
170  namelist /FILES/ ddkfile_1, ddkfile_2, ddkfile_3, wfkfile
171  namelist /PARAMETERS/ broadening, domega, maxomega, scissor, tolerance, do_antiresonant, do_temperature, &
172                        autoparal, max_ncpus
173  namelist /COMPUTATIONS/ num_lin_comp, lin_comp, num_nonlin_comp, nonlin_comp, &
174 &        num_linel_comp, linel_comp, num_nonlin2_comp, nonlin2_comp
175  namelist /TEMPERATURE/ epfile
176
177 ! *********************************************************************************
178
179 !Change communicator for I/O (mandatory!)
180  call abi_io_redirect(new_io_comm=xmpi_world)
181
182  call xmpi_init()
183  comm = xmpi_world
184  nproc = xmpi_comm_size(comm); my_rank = xmpi_comm_rank(comm)
185
186 !Initialize memory profiling if it is activated
187 !if a full abimem.mocc report is desired, set the argument of abimem_init to "2" instead of "0"
188 !note that abimem.mocc files can easily be multiple GB in size so don't use this option normally
189 #ifdef HAVE_MEM_PROFILING
190  call abimem_init(0)
191 #endif
192
193  call timein(tcpui,twall)
194  call timein(tcpui,twalli)
195  start_datetime = asctime()
196
197  if (my_rank == master) then
198    codename='OPTIC '//repeat(' ',18)
199    call herald(codename,abinit_version,std_out)
200
202    write(std_out,'(a)')' Please, give the name of the data file ...'
204    write(std_out,'(a)')' The name of the data file is :',filnam
205    write(std_out,'(a)')' Please, give the name of the output file ...'
207    write(std_out,'(a)')' The name of the output file is :',filnam_out
208    write(std_out,'(a)')' Please, give the root name for the (non)linear optical data output file ...'
210    write(std_out,'(a)')' The root name of the output files is :',trim(prefix)
211
213    if (open_file(filnam,msg,newunit=finunt,form='formatted') /= 0) then
214      MSG_ERROR(msg)
215    end if
216
217    write(msg,'(3a)') "From version 7.11.4, optic uses namelists as input.",ch10,&
218 &   "See e.g. ~/tests/tutorespfn/Input/toptic_2.in"
219    MSG_COMMENT(msg)
220
221    ! Setup some default values:
222    broadening = 1e-3_dp ! Ha
223    domega = 1e-3_dp ! Ha
224    maxomega = 1.0_dp ! Ha
225    scissor = 0.0_dp ! no scissor by default
226    tolerance = 1e-3_dp ! Ha
227    do_antiresonant = .TRUE. ! do use antiresonant approximation (only resonant transitions in the calculation)
228    do_temperature = .FALSE.
229
235    close(finunt)
236    ! Store filenames in array.
237    infiles = [wfkfile, ddkfile_1, ddkfile_2, ddkfile_3]
238
239    ! Validate input
240    if (num_nonlin_comp > 0 .and. all(nonlin_comp(1:num_nonlin_comp) == 0)) then
241      MSG_ERROR("nonlin_comp must be specified when num_nonlin_comp > 0")
242    end if
243    if (num_linel_comp > 0 .and. all(linel_comp(1:num_linel_comp) == 0)) then
244      MSG_ERROR("linel_comp must be specified when num_linel_comp > 0")
245    end if
246    if (num_nonlin2_comp > 0 .and. all(nonlin2_comp(1:num_nonlin2_comp) == 0)) then
247      MSG_ERROR("nonlin2_comp must be specified when num_nonlin2_comp > 0")
248    end if
249
250    ! Open the Wavefunction files
251    ! Note: Cannot use MPI-IO here because of prtwf=3.
252    ! If prtwf==3, the DDK file does not contain the wavefunctions but
253    ! this info is not reported in the header and the offsets in wfk_compute_offsets
254    ! are always computed assuming the presence of the cg
255
256    ! TODO: one should perform basic consistency tests for the GS WFK and the EVK files, e.g.
257    ! k-points and their order, spins, number of bands could differ in the four files.
258    ! Note indeed that we are assuming the same numer of bands in all the files.
259    call nctk_fort_or_ncfile(wfkfile, iomode, msg)
260    if (len_trim(msg) /= 0) MSG_ERROR(msg)
261    if (iomode == IO_MODE_MPI) iomode = IO_MODE_FORTRAN
263    ! Get header from the gs file
264    call hdr_copy(wfk0%hdr, hdr)
265
266    ! Read ddk here from WFK files or afterwards from EVK.nc
267    use_ncddk = .False.
268    do ii=1,3
269      use_ncddk(ii) = endswith(infiles(ii), "_EVK.nc")
270      if (.not. use_ncddk(ii)) then
271        call nctk_fort_or_ncfile(infiles(ii), iomode, msg)
272        if (len_trim(msg) /= 0) MSG_ERROR(msg)
273        if (iomode == IO_MODE_MPI) iomode = IO_MODE_FORTRAN
274        call wfk_open_read(wfks(ii), infiles(ii), formeig1, iomode, get_unit(), xmpi_comm_self)
275      end if
276    end do
277
278    ! Consistency check
279    do ii=1,2
280      if (.not. use_ncddk(ii) .and. .not. use_ncddk(ii+1)) then
281        if (wfk_compare(wfks(ii), wfks(ii+1)) /= 0) then
282          write(msg, "(2(a,i0,a))")"evkfile", ii," and evkfile ",ii+1, ", are not consistent. see above messages"
283          MSG_ERROR(msg)
284        end if
285      end if
286    end do
287
288    ep_nc_fname = 'test_EP.nc'; if (do_temperature) ep_nc_fname = epfile
289    do_ep_renorm = file_exists(ep_nc_fname)
290    ep_ntemp = 1
291    if (do_ep_renorm) then
292      call eprenorms_from_epnc(Epren,ep_nc_fname)
293      ep_ntemp = Epren%ntemp
294    else if (do_temperature) then
295      MSG_ERROR("You have asked for temperature but the epfile is not present !")
296    end if
297
298    ! autoparal section
299    if (autoparal /= 0 .and. max_ncpus /= 0) then
300      write(std_out,'(a)')"--- !Autoparal"
301      write(std_out,"(a)")'#Autoparal section for Optic runs.'
302      write(std_out,"(a)")   "info:"
303      write(std_out,"(a,i0)")"    autoparal: ",autoparal
304      write(std_out,"(a,i0)")"    max_ncpus: ",max_ncpus
305      write(std_out,"(a,i0)")"    nspinor: ",hdr%nspinor
306      write(std_out,"(a,i0)")"    nsppol: ",hdr%nsppol
307      write(std_out,"(a,i0)")"    nkpt: ",hdr%nkpt
308      write(std_out,"(a,i0)")"    mband: ",maxval(hdr%nband)
309
310      work_size = hdr%nkpt !* hdr%nsppol
311
312      ! List of configurations.
313      write(std_out,"(a)")"configurations:"
314      do ii=1,max_ncpus
315        if (ii > work_size) cycle
316        nks_per_proc = work_size / ii
317        nks_per_proc = nks_per_proc + MOD(work_size, ii)
318        eff = (one * work_size) / (ii * nks_per_proc)
319        write(std_out,"(a,i0)")"    - tot_ncpus: ",ii !* omp_ncpus
320        write(std_out,"(a,i0)")"      mpi_ncpus: ",ii
321        write(std_out,"(a,i0)")"      omp_ncpus: ",1
322        write(std_out,"(a,f12.9)")"      efficiency: ",eff
323        !write(,"(a,f12.2)")"      mem_per_cpu: ",mempercpu_mb
324      end do
325
326      write(std_out,'(a)')"..."
327      MSG_ERROR_NODUMP("aborting now")
328    end if
329
330  end if
331
332  call hdr_bcast(hdr,master,my_rank,comm)
333  !TODO put parameters in datastructure
335  call xmpi_bcast(domega,master,comm,ierr)
336  call xmpi_bcast(maxomega,master,comm,ierr)
337  call xmpi_bcast(scissor,master,comm,ierr)
338  call xmpi_bcast(tolerance,master,comm,ierr)
339  call xmpi_bcast(num_lin_comp,master,comm,ierr)
340  call xmpi_bcast(lin_comp,master,comm,ierr)
341  call xmpi_bcast(num_nonlin_comp,master,comm,ierr)
342  call xmpi_bcast(nonlin_comp,master,comm,ierr)
343  call xmpi_bcast(num_linel_comp,master,comm,ierr)
344  call xmpi_bcast(linel_comp,master,comm,ierr)
345  call xmpi_bcast(num_nonlin2_comp,master,comm,ierr)
346  call xmpi_bcast(nonlin2_comp,master,comm,ierr)
347  call xmpi_bcast(do_antiresonant,master,comm,ierr)
348  call xmpi_bcast(do_ep_renorm,master,comm,ierr)
349  call xmpi_bcast(ep_ntemp,master,comm,ierr)
350  if (do_ep_renorm) call eprenorms_bcast(Epren, master, comm)
351
352 !Extract info from the header
354  bantot=hdr%bantot
355  ecut=hdr%ecut_eff
356  natom=hdr%natom
357  nkpt=hdr%nkpt
358  nspinor=hdr%nspinor
359  nsppol=hdr%nsppol
360  ntypat=hdr%ntypat
361  occopt=hdr%occopt
362  rprimd(:,:)=hdr%rprimd(:,:)
363  ABI_ALLOCATE(nband,(nkpt*nsppol))
364  ABI_ALLOCATE(occ,(bantot))
365  !fermie=hdr%fermie
366  ! YG Fermi energy contained in the header of a NSCF computation is always 0 !!
367  occ(1:bantot)=hdr%occ(1:bantot)
368  nband(1:nkpt*nsppol)=hdr%nband(1:nkpt*nsppol)
369
370  nsym=hdr%nsym
371  ABI_ALLOCATE(symrel,(3,3,nsym))
372  ABI_ALLOCATE(symrec,(3,3,nsym))
373  symrel(:,:,:) = hdr%symrel(:,:,:)
374  do isym=1,nsym
375    call mati3inv(symrel(:,:,isym),symrec(:,:,isym))
376  end do
377
378  ABI_ALLOCATE(kpt,(3,nkpt))
379  kpt(:,:) = hdr%kptns(:,:)
380
381 !Get mband, as the maximum value of nband(nkpt)
382  mband=maxval(nband(:))
383  do ii=1,nkpt
384    if (nband(ii) /= mband) then
385      MSG_ERROR("nband must be constant across kpts")
386    end if
387  end do
388
389  ! Initializes crystal
390  call crystal_init(hdr%amu, cryst, 0, hdr%natom, hdr%npsp, hdr%ntypat, &
391 & hdr%nsym, rprimd, hdr%typat, hdr%xred, hdr%zionpsp, hdr%znuclpsp, 1, &
392 & (hdr%nspden==2 .and. hdr%nsppol==1),remove_inv, hdr%title,&
393 & symrel, hdr%tnons, hdr%symafm)
394
395  if (my_rank == master) then
396    write(std_out,*)
397    write(std_out,'(a,3f10.5,a)' )' rprimd(bohr)      =',rprimd(1:3,1)
398    write(std_out,'(a,3f10.5,a)' )'                    ',rprimd(1:3,2)
399    write(std_out,'(a,3f10.5,a)' )'                    ',rprimd(1:3,3)
400    write(std_out,'(a,i8)')       ' natom             =',natom
401    write(std_out,'(a,2i8)')      ' nkpt,mband        =',nkpt,mband
402    write(std_out,'(a, f10.5,a)' ) ' ecut              =',ecut,' Ha'
403  end if
404
405  ! Read the eigenvalues of ground-state and ddk files
406  ABI_ALLOCATE(eigen0,(mband*nkpt*nsppol))
407  ! MG: Do not understand why not [...,3]
408  ABI_ALLOCATE(eigen11,(2*mband*mband*nkpt*nsppol))
409  ABI_ALLOCATE(eigen12,(2*mband*mband*nkpt*nsppol))
410  ABI_ALLOCATE(eigen13,(2*mband*mband*nkpt*nsppol))
411
412  if (my_rank == master) then
413    ABI_ALLOCATE(eigtmp,(2*mband*mband))
414    ABI_ALLOCATE(eig0tmp,(mband))
415
416    do ii=1,3
417      if (.not. use_ncddk(ii)) cycle
418 #ifdef HAVE_NETCDF
420      varid = nctk_idname(ncid, "h1_matrix_elements")
421      outeig => eigen11
422      if (ii == 2) outeig => eigen12
423      if (ii == 3) outeig => eigen13
424      NCF_CHECK(nf90_get_var(ncid, varid, outeig, count=[2, mband, mband, nkpt, nsppol]))
425      NCF_CHECK(nf90_close(ncid))
426 #else
427      MSG_ERROR("Netcdf not available!")
428 #endif
429    end do
430
431    bdtot0_index=0 ; bdtot_index=0
432    do isppol=1,nsppol
433      do ikpt=1,nkpt
434        nband1=nband(ikpt+(isppol-1)*nkpt)
435        eigtmp = zero
436        eig0tmp = zero
437
439        eigen0(1+bdtot0_index:nband1+bdtot0_index)=eig0tmp(1:nband1)
440
441        ! Read DDK matrix elements from WFK
442        do ii=1,3
443          if (.not. use_ncddk(ii)) then
444            call wfk_read_eigk(wfks(ii), ikpt, isppol, xmpio_single, eigtmp)
445            if (ii == 1) eigen11(1+bdtot_index:2*nband1**2+bdtot_index)=eigtmp(1:2*nband1**2)
446            if (ii == 2) eigen12(1+bdtot_index:2*nband1**2+bdtot_index)=eigtmp(1:2*nband1**2)
447            if (ii == 3) eigen13(1+bdtot_index:2*nband1**2+bdtot_index)=eigtmp(1:2*nband1**2)
448            !ABI_CHECK(wfks(ii)%nband(ikpt,isppol) == nband1, "ddk1 nband1")
449          end if
450        end do
451        bdtot0_index=bdtot0_index+nband1
452        bdtot_index=bdtot_index+2*nband1**2
453      end do
454    end do
455
456    call wfk_close(wfk0)
457    do ii=1,3
458      if (.not. use_ncddk(ii)) call wfk_close(wfks(ii))
459    end do
460
461    ABI_DEALLOCATE(eigtmp)
462    ABI_DEALLOCATE(eig0tmp)
463  end if ! master
464
465  call xmpi_bcast(eigen0,master,comm,ierr)
466  call xmpi_bcast(eigen11,master,comm,ierr)
467  call xmpi_bcast(eigen12,master,comm,ierr)
468  call xmpi_bcast(eigen13,master,comm,ierr)
469
470 !---------------------------------------------------------------------------------
471 !gmet inversion
472  call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
473  call matr3inv(gmet,gmet_inv)
474
475 !---------------------------------------------------------------------------------
476 !derivative of occupation wrt the energy.
477  ABI_ALLOCATE(wtk,(nkpt))
478  wtk = hdr%wtk
479
480  ABI_ALLOCATE(doccde,(mband*nkpt*nsppol))
481
483  !WARNING no garantie that it works for other materials than insulators
484  nelect = hdr%nelect
485  tphysel = zero
486  ABI_ALLOCATE(istwfk,(nkpt))
487  ABI_ALLOCATE(npwarr,(nkpt))
488  istwfk = hdr%istwfk
489  npwarr = hdr%npwarr
490
491  call ebands_init(bantot, ks_ebands, nelect, doccde, eigen0, istwfk, kpt, &
492 & nband, nkpt, npwarr, nsppol, nspinor, tphysel, broadening, occopt, occ, wtk, &
493 & hdr%charge, hdr%kptopt, hdr%kptrlatt_orig, hdr%nshiftk_orig, hdr%shiftk_orig, &
494 & hdr%kptrlatt, hdr%nshiftk, hdr%shiftk)
495
496  !YG : should we use broadening for ebands_init
497  call ebands_update_occ(ks_ebands, -99.99d0)
498  fermie = ks_ebands%fermie
499  ABI_DEALLOCATE(istwfk)
500  ABI_DEALLOCATE(npwarr)
501
502 !---------------------------------------------------------------------------------
503 !size of the frequency range
504  nomega=int((maxomega+domega*0.001_dp)/domega)
505  maxomega = dble(nomega)*domega
506  ABI_ALLOCATE(cond_nd,(nomega))
507  ABI_ALLOCATE(cond_kg,(nomega))
508
509  optic_ncid = nctk_noid
510  if (my_rank == master) then
511    write(std_out,'(a,f10.5,a,f10.5,a)' )' fermie            =',fermie,' Ha',fermie*Ha_eV,' eV'
512    write(std_out,'(a,f10.5,a)')' Scissor shift     =', scissor, ' Ha'
513    write(std_out,'(a,f10.5,a)')' Tolerance on closeness to singularities     =', tolerance, ' Ha'
514    write(std_out,'(a,f10.5,a)')' Smearing factor      =', broadening, ' Ha'
515    if (do_antiresonant) then
516      write(std_out,'(a)') ' Will use the antiresonant approximation '
517    else
518      write(std_out,'(a)') ' Will not use the antiresonant approximation (only available for nonlin2 and linel components!) '
519    end if
520    write(std_out,'(a)') ' linear coeffs to be calculated : '
521    write(std_out,'(9i3)') lin_comp(1:num_lin_comp)
522    write(std_out,'(a)') ' non-linear coeffs to be calculated : '
523    write(std_out,'(27i4)') nonlin_comp(1:num_nonlin_comp)
524    write(std_out,'(a)') ' electronic part of electro-optic coeffs to be calculated :'
525    write(std_out,'(27i4)') linel_comp(1:num_linel_comp)
526    write(std_out,'(a)') ' non-linear coeffs (V2) to be calculated :'
527    write(std_out,'(27i4)') nonlin2_comp(1:num_nonlin2_comp)
528
529 #ifdef HAVE_NETCDF
530    ! Open netcdf file that will contain output results (only master is supposed to write)
531    NCF_CHECK_MSG(nctk_open_create(optic_ncid, strcat(prefix, "_OPTIC.nc"), xmpi_comm_self), "Creating _OPTIC.nc")
532
534    ! Note that we write the KS bands without EPH interaction (if any).
535    NCF_CHECK(hdr_ncwrite(hdr, optic_ncid, 666, nc_define=.True.))
536    NCF_CHECK(crystal_ncwrite(cryst, optic_ncid))
537    NCF_CHECK(ebands_ncwrite(ks_ebands, optic_ncid))
538
539    ! Add optic input variables.
540    NCF_CHECK(nctk_def_dims(optic_ncid, [nctkdim_t("ntemp", ep_ntemp), nctkdim_t("nomega", nomega)], defmode=.True.))
541
542    ncerr = nctk_def_iscalars(optic_ncid, [character(len=nctk_slen) :: "do_antiresonant", "do_ep_renorm"])
543    NCF_CHECK(ncerr)
544    ncerr = nctk_def_dpscalars(optic_ncid, [character(len=nctk_slen) :: &
545    "broadening", "domega", "maxomega", "scissor", "tolerance"])
546    NCF_CHECK(ncerr)
547
548    ! Define arrays containing output results
549    ncerr = nctk_def_arrays(optic_ncid, [ &
550    nctkarr_t('wmesh', "dp", "nomega") &
551    ])
552    NCF_CHECK(ncerr)
553
554    if (num_lin_comp > 0) then
555      ! Linear optic results.
556      NCF_CHECK(nctk_def_dims(optic_ncid, nctkdim_t("linopt_ncomp", num_lin_comp)))
557      ncerr = nctk_def_arrays(optic_ncid, [ &
558      nctkarr_t('linopt_components', "int", "linopt_ncomp"), &
559      nctkarr_t('linopt_epsilon', "dp", "two, nomega, linopt_ncomp, ntemp") &
560      ])
561      NCF_CHECK(ncerr)
562    end if
563
564    if (num_nonlin_comp > 0) then
565      ! Second harmonic generation.
566      NCF_CHECK(nctk_def_dims(optic_ncid, nctkdim_t("shg_ncomp", num_nonlin_comp)))
567      ncerr = nctk_def_arrays(optic_ncid, [ &
568      nctkarr_t('shg_components', "int", "shg_ncomp"), &
569      nctkarr_t('shg_inter2w', "dp", "two, nomega, shg_ncomp, ntemp"), &
570      nctkarr_t('shg_inter1w', "dp", "two, nomega, shg_ncomp, ntemp"), &
571      nctkarr_t('shg_intra2w', "dp", "two, nomega, shg_ncomp, ntemp"), &
572      nctkarr_t('shg_intra1w', "dp", "two, nomega, shg_ncomp, ntemp"), &
573      nctkarr_t('shg_intra1wS', "dp", "two, nomega, shg_ncomp, ntemp"), &
574      nctkarr_t('shg_chi2tot', "dp", "two, nomega, shg_ncomp, ntemp") &
575      ])
576      NCF_CHECK(ncerr)
577    end if
578
579    if (num_linel_comp > 0) then
580      ! linear electro-optic (LEO) susceptibility
581      NCF_CHECK(nctk_def_dims(optic_ncid, nctkdim_t("leo_ncomp", num_linel_comp)))
582      ncerr = nctk_def_arrays(optic_ncid, [ &
583      nctkarr_t('leo_components', "int", "leo_ncomp"), &
584      nctkarr_t('leo_chi', "dp", "two, nomega, leo_ncomp, ntemp"), &
585      nctkarr_t('leo_eta', "dp", "two, nomega, leo_ncomp, ntemp"), &
586      nctkarr_t('leo_sigma', "dp", "two, nomega, leo_ncomp, ntemp"), &
587      nctkarr_t('leo_chi2tot', "dp", "two, nomega, leo_ncomp, ntemp") &
588      ])
589      NCF_CHECK(ncerr)
590    end if
591
592    if (num_nonlin2_comp > 0) then
593      ! non-linear electro-optic susceptibility
594      NCF_CHECK(nctk_def_dims(optic_ncid, nctkdim_t("leo2_ncomp", num_nonlin2_comp)))
595      ncerr = nctk_def_arrays(optic_ncid, [ &
596      nctkarr_t('leo2_components', "int", "leo2_ncomp"), &
597      nctkarr_t('leo2_chiw', "dp", "two, nomega, leo2_ncomp, ntemp"), &
598      nctkarr_t('leo2_etaw', "dp", "two, nomega, leo2_ncomp, ntemp"), &
599      nctkarr_t('leo2_chi2w', "dp", "two, nomega, leo2_ncomp, ntemp"), &
600      nctkarr_t('leo2_eta2w', "dp", "two, nomega, leo2_ncomp, ntemp"), &
601      nctkarr_t('leo2_sigmaw', "dp", "two, nomega, leo2_ncomp, ntemp"), &
602      nctkarr_t('leo2_chi2tot', "dp", "two, nomega, leo2_ncomp, ntemp") &
603      ])
604      NCF_CHECK(ncerr)
605    end if
606
607    NCF_CHECK(nctk_set_datamode(optic_ncid))
608
609    ! Write wmesh here.
610    ABI_ALLOCATE(wmesh, (nomega))
611    do ii=1,nomega
612      ! This to be consistent with the value used in m_optic_tools
613      ! In principle wmesh should be passed to the children and a lot of code
614      ! should be rewritten to be more cache-friendly ...
615      wmesh(ii) = (ii-1)*domega * (13.60569172*2._dp)
616    end do
617    NCF_CHECK(nf90_put_var(optic_ncid, nctk_idname(optic_ncid, "wmesh"), wmesh))
618    ABI_FREE(wmesh)
619
620    if (num_lin_comp > 0) then
621      NCF_CHECK(nf90_put_var(optic_ncid, nctk_idname(optic_ncid, "linopt_components"), lin_comp(1:num_lin_comp)))
622    end if
623    if (num_nonlin_comp > 0) then
624      NCF_CHECK(nf90_put_var(optic_ncid, nctk_idname(optic_ncid, "shg_components"), nonlin_comp(1:num_nonlin_comp)))
625    end if
626    if (num_linel_comp > 0) then
627      NCF_CHECK(nf90_put_var(optic_ncid, nctk_idname(optic_ncid, "leo_components"), linel_comp(1:num_linel_comp)))
628    end if
629    if (num_nonlin2_comp > 0) then
630      NCF_CHECK(nf90_put_var(optic_ncid, nctk_idname(optic_ncid, "leo2_components"), nonlin2_comp(1:num_nonlin2_comp)))
631    end if
632
633    ! Write optic input variables.
634    ii = 0; if (do_antiresonant) ii = 1
635    jj = 0; if (do_ep_renorm) jj = 1
636    ncerr = nctk_write_iscalars(optic_ncid, [character(len=nctk_slen) :: &
637    "do_antiresonant", "do_ep_renorm"], &
638    [ii, jj])
639    NCF_CHECK(ncerr)
640
641    ncerr = nctk_write_dpscalars(optic_ncid, [character(len=nctk_slen) :: &
642    "broadening", "domega", "maxomega", "scissor", "tolerance"], &
643    [broadening, domega, maxomega, scissor, tolerance])
644    NCF_CHECK(ncerr)
645 #endif
646  end if
647
648  ABI_ALLOCATE(symcart,(3,3,nsym))
649  !YG : we need to transpose gprimd since matrinv give the transpose of the inverse !
650  gprimd_trans = transpose(gprimd)
651  call sym2cart(gprimd_trans,nsym,rprimd,symrel,symcart)
652
653  ABI_ALLOCATE(pmat,(mband,mband,nkpt,3,nsppol))
654  call wrtout(std_out," optic : Call pmat2cart","COLL")
655
656  call pmat2cart(eigen11,eigen12,eigen13,mband,nkpt,nsppol,pmat,rprimd)
657  call pmat_renorm(fermie, eigen0, mband, nkpt, nsppol, pmat, scissor)
658
659 !IN CALLED ROUTINE
660 !call linopt(nspin,,nkpt,wkpt,nsymcrys,symcrys,nstval,occv,evalv,efermi,pmat,v1,v2,nmesh,de,scissor,brod)
661 !
662 !v1,v2=desired component of the dielectric function(integer) 1=x,2=y,3=z
663 !nmesh=desired number of energy mesh points(integer)
664 !de=desired step in energy(real); nmesh*de=maximum energy
665 !scissor=scissors shift in Ha(real)
667
668  ! optical frequency dependent dielectric function for semiconductors
669  call wrtout(std_out," optic : Call linopt","COLL")
670
671  do itemp=1,ep_ntemp
672    call ebands_copy(ks_ebands, eph_ebands)
673    if (do_ep_renorm) call renorm_bst(Epren, eph_ebands, cryst, itemp, do_lifetime=.True.,do_check=.True.)
674    do ii=1,num_lin_comp
675      lin1 = int(lin_comp(ii)/10.0_dp)
676      lin2 = mod(lin_comp(ii),10)
677      write(msg,*) ' linopt ', lin1,lin2
678      call wrtout(std_out,msg,"COLL")
679      call int2char4(lin1,s1)
680      call int2char4(lin2,s2)
681      call int2char4(itemp,stemp)
682      ABI_CHECK((s1(1:1)/='#'),'Bug: string length too short!')
683      ABI_CHECK((s2(1:1)/='#'),'Bug: string length too short!')
684      ABI_CHECK((stemp(1:1)/='#'),'Bug: string length too short!')
686      if (do_ep_renorm) tmp_radix = trim(prefix)//"_"//trim(s1)//"_"//trim(s2)//"_T"//trim(stemp)
687      call linopt(ii,itemp,nsppol,ucvol,nkpt,wtk,nsym,symcart,mband,ks_ebands,eph_ebands,fermie,pmat, &
689    end do
690    call ebands_free(eph_ebands)
691  end do
692
693  if (do_ep_renorm) call eprenorms_free(Epren)
694
695  ! second harmonic generation susceptibility for semiconductors
696  call wrtout(std_out," optic : Call nlinopt","COLL")
697  do ii=1,num_nonlin_comp
698    nlin1 = int( nonlin_comp(ii)/100.0_dp)
699    nlin2 = int((nonlin_comp(ii)-nlin1*100.0_dp)/10.0_dp)
700    nlin3 = mod( nonlin_comp(ii),10)
701    write(msg,*) ' nlinopt ', nlin1,nlin2,nlin3
702    call wrtout(std_out,msg,"COLL")
703    call int2char4(nlin1,s1)
704    call int2char4(nlin2,s2)
705    call int2char4(nlin3,s3)
706    ABI_CHECK((s1(1:1)/='#'),'Bug: string length too short!')
707    ABI_CHECK((s2(1:1)/='#'),'Bug: string length too short!')
708    ABI_CHECK((s3(1:1)/='#'),'Bug: string length too short!')
710    itemp = 1
711    call nlinopt(ii,itemp,nsppol,ucvol,nkpt,wtk,nsym,symcart,mband,eigen0,fermie,pmat,&
713  end do
714
715  ! linear electro-optic susceptibility for semiconductors
716  call wrtout(std_out," optic : Call linelop","COLL")
717  do ii=1,num_linel_comp
718    linel1 = int( linel_comp(ii)/100.0_dp)
719    linel2 = int((linel_comp(ii)-linel1*100.0_dp)/10.0_dp)
720    linel3 = mod( linel_comp(ii),10)
721    write(msg,*) ' linelop ',linel1,linel2,linel3
722    call wrtout(std_out,msg,"COLL")
723    call int2char4(linel1,s1)
724    call int2char4(linel2,s2)
725    call int2char4(linel3,s3)
727    itemp = 1
728    call linelop(ii,itemp,nsppol,ucvol,nkpt,wtk,nsym,symcart,mband,eigen0,occ,fermie,pmat,&
730  end do
731
732  call wrtout(std_out," optic : Call nonlinopt","COLL")
733  do ii=1,num_nonlin2_comp
734    nonlin1 = int( nonlin2_comp(ii)/100.0_dp)
735    nonlin2 = int((nonlin2_comp(ii)-nonlin1*100.0_dp)/10.0_dp)
736    nonlin3 = mod( nonlin2_comp(ii),10)
737    write(msg,*) ' nonlinopt ',nonlin1,nonlin2,nonlin3
738    call wrtout(std_out,msg,"COLL")
739    call int2char4(nonlin1,s1)
740    call int2char4(nonlin2,s2)
741    call int2char4(nonlin3,s3)
743    itemp = 1
744    call nonlinopt(ii,itemp,nsppol,ucvol,nkpt,wtk,nsym,symcart,mband,eigen0,occ,fermie,pmat,&
746  end do
747
748  ABI_DEALLOCATE(nband)
749  ABI_DEALLOCATE(occ)
750  ABI_DEALLOCATE(eigen11)
751  ABI_DEALLOCATE(eigen12)
752  ABI_DEALLOCATE(eigen13)
753  ABI_DEALLOCATE(eigen0)
754  ABI_DEALLOCATE(doccde)
755  ABI_DEALLOCATE(wtk)
756  ABI_DEALLOCATE(cond_nd)
757  ABI_DEALLOCATE(cond_kg)
758  ABI_DEALLOCATE(kpt)
759  ABI_DEALLOCATE(symrel)
760  ABI_DEALLOCATE(symrec)
761  ABI_DEALLOCATE(symcart)
762  ABI_DEALLOCATE(pmat)
763
764  call hdr_free(hdr)
765  call ebands_free(ks_ebands)
766  call crystal_free(cryst)
767
768  call timein(tcpu,twall)
769
770  tsec(1)=tcpu-tcpui
771  tsec(2)=twall-twalli
772
773  if (my_rank == master) then
774    write(std_out,'(a,80a,a,a,a)' )ch10,('=',ii=1,80),ch10,ch10,' Calculation completed.'
775    write(std_out, '(a,a,a,f13.1,a,f13.1)' ) &
776 &   '-',ch10,'- Proc.   0 individual time (sec): cpu=',tsec(1),'  wall=',tsec(2)
777  end if
778
779  call xmpi_sum(tsec,comm,ierr)
780
781  if (my_rank == master) then
782    ! Write YAML document with the final summary.
783    ! we use this doc to test whether the calculation is completed.
784    write(std_out,"(a)")"--- !FinalSummary"
785    write(std_out,"(a)")"program: optic"
786    write(std_out,"(2a)")"version: ",trim(abinit_version)
787    write(std_out,"(2a)")"start_datetime: ",start_datetime
788    write(std_out,"(2a)")"end_datetime: ",asctime()
789    write(std_out,"(a,f13.1)")"overall_cpu_time: ",tsec(1)
790    write(std_out,"(a,f13.1)")"overall_wall_time: ",tsec(2)
791    write(std_out,"(a,i0)")"mpi_procs: ",xmpi_comm_size(xmpi_world)
793    !write(std_out,"(a,i0)")"num_warnings: ",nwarning
795    write(std_out,"(a)")"..."
796    call flush_unit(std_out)
797  end if
798
799 #ifdef HAVE_NETCDF
800  if (my_rank == master) then
801    NCF_CHECK(nf90_close(optic_ncid))
802  end if
803 #endif
804
805 !Write information on file about the memory before ending mpi module, if memory profiling is enabled
806  call abinit_doctor(filnam)
807
808  call xmpi_end()
809
810  end program optic
```