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