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), ABI_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 ! Ha, 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 : !!***
|
