Line data Source code
1 : !!****m* ABINIT/m_rttddft_output
2 : !! NAME
3 : !! m_rttddft_ouptut
4 : !!
5 : !! FUNCTION
6 : !! Manages most output of RT-TDDFT runs
7 : !!
8 : !! COPYRIGHT
9 : !! Copyright (C) 2021-2026 ABINIT group (FB)
10 : !! This file is distributed under the terms of the
11 : !! GNU General Public License, see ~abinit/COPYING
12 : !! or http://www.gnu.org/copyleft/gpl.txt .
13 : !!
14 : !! SOURCE
15 :
16 : #if defined HAVE_CONFIG_H
17 : #include "config.h"
18 : #endif
19 :
20 : #include "abi_common.h"
21 :
22 : module m_rttddft_output
23 :
24 : use netcdf
25 :
26 : use defs_basis
27 : use defs_abitypes, only: MPI_type
28 : use defs_datatypes, only: pseudopotential_type
29 : use m_common, only: prteigrs
30 : use m_crystal, only: crystal_t
31 : use m_dtfil, only: datafiles_type
32 : use m_dtset, only: dataset_type
33 : use m_ebands, only: ebands_t
34 : use m_epjdos, only: epjdos_t
35 : use m_errors, only: msg_hndl, assert
36 : use m_ioarr, only: fftdatar_write
37 : use m_io_tools, only: open_file, close_unit
38 : use m_iowf, only: outwf
39 : use m_mpinfo, only: iwrite_fftdatar
40 : use m_paral_atom, only: get_my_atmtab, free_my_atmtab
41 : use m_profiling_abi, only: abimem_record
42 : use m_rttddft_tdks, only: tdks_type
43 : use m_specialmsg, only: wrtout
44 : use m_xmpi, only: xmpi_comm_rank
45 :
46 : implicit none
47 :
48 : private
49 : !!***
50 :
51 : public :: rttddft_output
52 : !!***
53 :
54 : contains
55 : !!***
56 :
57 : !!****f* m_rttddft_output/rttddft_output
58 : !!
59 : !! NAME
60 : !! rttddft_output
61 : !!
62 : !! FUNCTION
63 : !! Main output subroutine
64 : !!
65 : !! INPUTS
66 : !! dtfil <type datafiles_type> = infos about file names, file unit numbers
67 : !! dtset <type(dataset_type)> = all input variables for this dataset
68 : !! istep <integer> = step number
69 : !! mpi_enreg <MPI_type> = MPI-parallelisation information
70 : !! psps <type(pseudopotential_type)> = variables related to pseudopotentials
71 : !! tdks <type(tdks_type)> = the tdks object to initialize
72 : !!
73 : !! OUTPUT
74 : !!
75 : !! SOURCE
76 1341 : subroutine rttddft_output(dtfil, dtset, istep, mpi_enreg, psps, tdks)
77 :
78 : !Arguments ------------------------------------
79 : !scalars
80 : integer, intent(in) :: istep
81 : type(datafiles_type), intent(inout) :: dtfil
82 : type(dataset_type), intent(inout) :: dtset
83 : type(MPI_type), intent(inout) :: mpi_enreg
84 : type(pseudopotential_type), intent(inout) :: psps
85 : type(tdks_type), intent(inout) :: tdks
86 : !arrays
87 :
88 : !Local variables-------------------------------
89 : !scalars
90 : character(len=500) :: msg
91 : character(len=100) :: fmt
92 : character(len=20) :: access
93 : integer :: i
94 : !arrays
95 :
96 : ! *************************************************************************
97 :
98 : !** Special case of first step
99 1341 : if (istep == tdks%first_step) then
100 46 : access = "sequential"
101 46 : if (dtset%td_restart>0) access = "append"
102 :
103 : ! Open energy file and writes header if needed
104 46 : if (mpi_enreg%me == 0) then
105 18 : if (open_file(tdks%fname_tdener,msg,newunit=tdks%tdener_unit,status='unknown',form='formatted',access=access) /= 0) ABI_ERROR(msg)
106 : end if
107 46 : if (dtset%td_restart==0) then
108 26 : write(msg,'(a)') "# RT-TDDFT -- Energy file. All quantities are in Hartree atomic units."
109 26 : call wrtout(tdks%tdener_unit,msg)
110 : write(msg,'(a)') "# step time E_total E_kinetic E_hartree E_xc E_ewald &
111 26 : & E_corepsp E_localpsp E_nonlocalpsp E_paw E_entropy E_vdw"
112 26 : call wrtout(tdks%tdener_unit,msg)
113 : end if
114 :
115 : ! Open electric field file and writes header if needed
116 46 : if (dtset%td_ef_type /= 0) then
117 21 : if (mpi_enreg%me == 0) then
118 7 : if (open_file(tdks%fname_tdef,msg,newunit=tdks%tdef_unit,status='unknown',form='formatted',access=access) /= 0) ABI_ERROR(msg)
119 : end if
120 21 : if (dtset%td_restart==0) then
121 11 : write(msg,'(a)') "# RT-TDDFT -- Electric field file. All quantities are in Hartree atomic units."
122 11 : call wrtout(tdks%tdef_unit,msg)
123 11 : write(msg,'(a)') "# step time E_x E_y E_z A_x A_y A_z A_ext_x A_ext_y A_ext_z A_ind_x A_ind_y A_ind_z"
124 11 : call wrtout(tdks%tdef_unit,msg)
125 : end if
126 : end if
127 :
128 : ! Open current file and writes header if needed
129 46 : if (dtset%prtcurrent /= 0) then
130 21 : if (mpi_enreg%me == 0) then
131 7 : if (open_file(tdks%fname_current,msg,newunit=tdks%current_unit,status='unknown',form='formatted',access=access) /= 0) ABI_ERROR(msg)
132 : end if
133 21 : if (dtset%td_restart==0) then
134 11 : write(msg,'(a)') "# RT-TDDFT -- Current density file. All quantities are in Hartree atomic units."
135 11 : call wrtout(tdks%current_unit,msg)
136 11 : write(msg,'(a)') "# step time J_x J_y J_z"
137 11 : call wrtout(tdks%current_unit,msg)
138 : end if
139 : end if
140 : end if
141 :
142 : !** Writes some info in main output file
143 1341 : write(msg,'(a,a,f14.6,a)') ch10,'Total energy = ', tdks%etot,' Ha'
144 1341 : call wrtout(ab_out,msg)
145 1341 : if (do_write_log) call wrtout(std_out,msg)
146 :
147 1341 : write(msg,'(a,f18.10,a)') 'Integrated density (ie. total nb of electrons) = ', &
148 21767625 : & sum(tdks%rhor(:,1))*tdks%ucvol/tdks%nfftf, ch10
149 1341 : call wrtout(ab_out,msg)
150 1341 : if (do_write_log) call wrtout(std_out,msg)
151 :
152 : !** Writes in energy file
153 1341 : write(msg,'(i0,1X,f15.5,11(f14.8,1X))') istep-1, (istep-1)*tdks%dt, tdks%etot, tdks%energies%e_kinetic, &
154 1341 : & tdks%energies%e_hartree, tdks%energies%e_xc, tdks%energies%e_ewald, &
155 1341 : & tdks%energies%e_corepsp, tdks%energies%e_localpsp, tdks%energies%e_nlpsp_vfock, &
156 2682 : & tdks%energies%paw%epaw, tdks%energies%e_entropy, tdks%energies%e_vdw_dftd
157 1341 : call wrtout(tdks%tdener_unit,msg)
158 :
159 : !** Writes TD elec. field and associated vector potential if needed
160 1341 : if (dtset%td_ef_type /= 0) then
161 1100 : write(msg,'(i0,1X,f15.5,1X,12(f14.8,1X))') istep, istep*tdks%dt, tdks%tdef%efield(:), tdks%tdef%vecpot(:), &
162 2200 : & tdks%tdef%vecpot_ext(:), tdks%tdef%vecpot_ind(:,1)
163 1100 : call wrtout(tdks%tdef_unit,msg)
164 : end if
165 :
166 : !** Writes TD current density if needed
167 1341 : if (dtset%prtcurrent /= 0) then
168 1100 : if (dtset%nsppol == 1) then
169 1100 : fmt = '(i0,1X,f15.5,1X,3(f14.8,1X))'
170 : else
171 0 : fmt = '(i0,1X,f15.5,1X,3(f14.8,1X),3(f14.8,1X))'
172 : end if
173 2200 : write(msg,fmt) istep, istep*tdks%dt, (tdks%current(:,i),i=1,dtset%nsppol)
174 1100 : call wrtout(tdks%current_unit,msg)
175 : end if
176 :
177 : !** Writes additional optional properties
178 : !Update header, with evolving variables
179 : call tdks%hdr%update(tdks%bantot,tdks%etot,tdks%energies%e_fermie,tdks%energies%e_fermih, &
180 : & tdks%hdr%residm,tdks%rprimd,tdks%occ0,tdks%pawrhoij, &
181 : & tdks%xred,dtset%amu_orig,comm_atom=mpi_enreg%comm_atom, &
182 1341 : & mpi_atmtab=mpi_enreg%my_atmtab)
183 :
184 : !Computed at actual step
185 1341 : if (mod(istep,dtset%td_prtstr) == 0) then
186 1341 : call prt_den(dtfil,dtset,istep,mpi_enreg,psps,tdks)
187 1341 : if (dtset%prtwf > 0) then
188 0 : call prt_wfk(dtfil,dtset,istep,mpi_enreg,psps,tdks)
189 0 : call prt_restart(dtfil,dtset,istep,mpi_enreg,tdks)
190 : end if
191 : end if
192 :
193 : !Computed at previous step
194 1341 : if (mod(istep-1,dtset%td_prtstr) == 0) then
195 1341 : call prt_eig(dtfil,dtset,istep-1,mpi_enreg,tdks)
196 1341 : call prt_occ(dtfil,dtset,istep-1,mpi_enreg,tdks)
197 1341 : call prt_dos(dtfil,dtset,istep-1,mpi_enreg,psps,tdks)
198 : end if
199 :
200 : !** Special case of last step
201 1341 : if (istep == tdks%first_step+tdks%ntime-1) then
202 46 : if (mod(istep,dtset%td_prtstr) /= 0 .or. dtset%prtwf <= 0) then
203 46 : call prt_wfk(dtfil,dtset,istep,mpi_enreg,psps,tdks,force_write=.TRUE.)
204 46 : call prt_restart(dtfil,dtset,istep,mpi_enreg,tdks)
205 : end if
206 46 : if (mpi_enreg%me == 0) then
207 : !close all files
208 18 : if (close_unit(tdks%tdener_unit,msg) /= 0) ABI_ERROR(msg)
209 18 : if (close_unit(tdks%tdrestart_unit,msg) /= 0) ABI_ERROR(msg)
210 18 : if (dtset%td_ef_type /= 0) then
211 7 : if (close_unit(tdks%tdef_unit,msg) /= 0) ABI_ERROR(msg)
212 : end if
213 18 : if (dtset%prtcurrent /= 0) then
214 7 : if (close_unit(tdks%current_unit,msg) /= 0) ABI_ERROR(msg)
215 : end if
216 : end if
217 : end if
218 :
219 1341 : end subroutine rttddft_output
220 : !!***
221 :
222 : !!****f* m_rttddft_output/prt_eig
223 : !!
224 : !! NAME
225 : !! prt_eig
226 : !!
227 : !! FUNCTION
228 : !! Outputs eigenvalues
229 : !!
230 : !! INPUTS
231 : !! dtfil <type datafiles_type> = infos about file names, file unit numbers
232 : !! dtset <type(dataset_type)> = all input variables for this dataset
233 : !! istep <integer> = step number
234 : !! mpi_enreg <MPI_type> = MPI-parallelisation information
235 : !! tdks <type(tdks_type)> = the tdks object to initialize
236 : !!
237 : !! OUTPUT
238 : !!
239 : !! SIDE EFFECTS
240 : !!
241 : !! SOURCE
242 1341 : subroutine prt_eig(dtfil, dtset, istep, mpi_enreg, tdks)
243 :
244 : !Arguments ------------------------------------
245 : !scalars
246 : integer, intent(in) :: istep
247 : type(datafiles_type), intent(inout) :: dtfil
248 : type(dataset_type), intent(inout) :: dtset
249 : type(MPI_type), intent(inout) :: mpi_enreg
250 : type(tdks_type), intent(inout) :: tdks
251 : !arrays
252 :
253 : !Local variables-------------------------------
254 : !scalars
255 : integer,parameter :: enunit=0, option=3
256 : integer :: me
257 : integer :: spacecomm
258 : real(dp) :: vxcavg_dum
259 : character(len=fnlen) :: fname
260 : character(len=24) :: step_nb
261 : !arrays
262 2682 : real(dp) :: resid(dtset%mband*dtset%nkpt*dtset%nsppol)
263 :
264 : ! *************************************************************************
265 :
266 1341 : spacecomm = mpi_enreg%comm_cell
267 1341 : me = xmpi_comm_rank(spacecomm)
268 :
269 : !to avoid some useless printing in the output file when both prteig and prtvol are set to 0
270 1341 : if (dtset%prteig /= 0 .or. dtset%prtvol /= 0) then
271 :
272 0 : write(step_nb,*) istep
273 0 : fname = trim(dtfil%filnam_ds(4))//'_'//trim(adjustl(step_nb))//'_EIG'
274 0 : resid = zero
275 0 : vxcavg_dum=zero
276 :
277 0 : if(me==0)then
278 : call prteigrs(tdks%eigen,enunit,tdks%energies%e_fermie,tdks%energies%e_fermih, &
279 : & fname,ab_out,dtset%iscf,dtset%kptns,dtset%kptopt,dtset%mband, &
280 : & dtset%nband,dtset%nbdbuf,dtset%nkpt,0,dtset%nsppol,tdks%occ0, &
281 : & dtset%occopt,option,dtset%prteig,dtset%prtvol,resid,dtset%tolwfr, &
282 0 : & vxcavg_dum,dtset%wtk)
283 : end if
284 : end if
285 :
286 1341 : end subroutine prt_eig
287 : !!***
288 :
289 : !!****f* m_rttddft_output/prt_occ
290 : !!
291 : !! NAME
292 : !! prt_occ
293 : !!
294 : !! FUNCTION
295 : !! Outputs occupation numbers
296 : !!
297 : !! INPUTS
298 : !! dtfil <type datafiles_type> = infos about file names, file unit numbers
299 : !! dtset <type(dataset_type)> = all input variables for this dataset
300 : !! istep <integer> = step number
301 : !! mpi_enreg <MPI_type> = MPI-parallelisation information
302 : !! tdks <type(tdks_type)> = the tdks object to initialize
303 : !!
304 : !! OUTPUT
305 : !!
306 : !! SOURCE
307 1341 : subroutine prt_occ(dtfil, dtset, istep, mpi_enreg, tdks)
308 :
309 : !Arguments ------------------------------------
310 : !scalars
311 : integer, intent(in) :: istep
312 : type(datafiles_type), intent(inout) :: dtfil
313 : type(dataset_type), intent(inout) :: dtset
314 : type(MPI_type), intent(inout) :: mpi_enreg
315 : type(tdks_type), intent(inout) :: tdks
316 : !arrays
317 :
318 : !Local variables-------------------------------
319 : !scalars
320 : integer :: band_index
321 : integer :: iband, ii, ikpt, isppol
322 : integer :: me
323 : integer :: nkpt
324 : integer :: nband_k, nsppol
325 : integer :: temp_unit
326 : !arrays
327 : character(len=fnlen) :: fname
328 : character(len=4) :: ibnd_fmt, ikpt_fmt
329 : character(len=500) :: msg
330 : character(len=24) :: step_nb
331 :
332 :
333 : ! *************************************************************************
334 :
335 1341 : if (dtset%prtocc > 0) then
336 0 : me = xmpi_comm_rank(mpi_enreg%comm_cell)
337 :
338 0 : write(step_nb,*) istep
339 0 : fname = trim(dtfil%filnam_ds(4))//'_'//trim(adjustl(step_nb))//'_OCC'
340 :
341 0 : if (open_file(fname, msg, newunit=temp_unit, status='unknown', form='formatted') /= 0) then
342 0 : ABI_ERROR(msg)
343 : end if
344 :
345 0 : nkpt = dtset%nkpt
346 0 : nsppol = dtset%nsppol
347 :
348 0 : if(me==0)then
349 : band_index=0
350 0 : do isppol=1,nsppol
351 :
352 0 : if(nsppol==2)then
353 0 : if(isppol==1)write(msg, '(2a)' ) ch10,' SPIN UP channel '
354 0 : if(isppol==2)write(msg, '(2a)' ) ch10,' SPIN DOWN channel '
355 0 : call wrtout(temp_unit,msg)
356 : end if
357 0 : ikpt_fmt="i4" ; if(nkpt>=10000)ikpt_fmt="i6" ; if(nkpt>=1000000)ikpt_fmt="i9"
358 0 : if (nsppol==2.and.isppol==1) then
359 : write(msg, '(a,'//ikpt_fmt//',2x,a)' ) &
360 0 : 'Occupation numbers for nkpt=',nkpt,'k points, SPIN UP:'
361 0 : else if (nsppol==2.and.isppol==2) then
362 : write(msg, '(a,'//ikpt_fmt//',2x,a)' ) &
363 0 : 'Occupation numbers for nkpt=',nkpt,'k points, SPIN DOWN:'
364 : else
365 : write(msg, '(a,'//ikpt_fmt//',2x,a)' ) &
366 0 : 'Occupation numbers for nkpt=',nkpt,'k points:'
367 : end if
368 0 : call wrtout(temp_unit,msg)
369 0 : do ikpt=1,nkpt
370 0 : nband_k=dtset%nband(ikpt+(isppol-1)*nkpt)
371 0 : ikpt_fmt="i4" ; if(nkpt>=10000)ikpt_fmt="i6" ; if(nkpt>=1000000)ikpt_fmt="i9"
372 0 : ibnd_fmt="i3" ; if(nband_k>=1000)ibnd_fmt="i6" ; if(nband_k>=1000000)ibnd_fmt="i9"
373 : write(msg, '(a,'//ikpt_fmt//',a,'//ibnd_fmt//',a,f9.5,a,3f8.4,a)' ) &
374 0 : & ' kpt#',ikpt,', nband=',nband_k,', wtk=',dtset%wtk(ikpt)+tol10,', kpt=',&
375 0 : & dtset%kptns(1:3,ikpt)+tol10,' (reduced coord)'
376 0 : call wrtout(temp_unit,msg)
377 0 : do ii=0,(nband_k-1)/6
378 0 : write(msg, '(1p,6e12.4)')(tdks%occ(iband+band_index),iband=1+6*ii,min(6+6*ii,nband_k))
379 0 : call wrtout(temp_unit,msg)
380 : end do
381 0 : band_index=band_index+nband_k
382 : end do
383 : end do
384 : end if
385 :
386 0 : close(temp_unit)
387 : end if
388 :
389 1341 : end subroutine prt_occ
390 : !!***
391 :
392 : !!****f* m_rttddft_output/prt_den
393 : !!
394 : !! NAME
395 : !! prt_den
396 : !!
397 : !! FUNCTION
398 : !! Outputs the electronic density
399 : !!
400 : !! INPUTS
401 : !! dtfil <type datafiles_type> = infos about file names, file unit numbers
402 : !! dtset <type(dataset_type)> = all input variables for this dataset
403 : !! istep <integer> = step number
404 : !! mpi_enreg <MPI_type> = MPI-parallelisation information
405 : !! psps <type(pseudopotential_type)> = variables related to pseudopotentials
406 : !! tdks <type(tdks_type)> = the tdks object to initialize
407 : !!
408 : !! OUTPUT
409 : !!
410 : !! SOURCE
411 1341 : subroutine prt_den(dtfil, dtset, istep, mpi_enreg, psps, tdks)
412 :
413 : !Arguments ------------------------------------
414 : !scalars
415 : integer, intent(in) :: istep
416 : type(datafiles_type), intent(inout) :: dtfil
417 : type(dataset_type), intent(inout) :: dtset
418 : type(MPI_type), intent(inout) :: mpi_enreg
419 : type(pseudopotential_type), intent(inout) :: psps
420 : type(tdks_type), intent(inout) :: tdks
421 : !arrays
422 :
423 : !Local variables-------------------------------
424 : !scalars
425 : integer,parameter :: cplex1=1
426 : integer :: bantot
427 : integer :: iatom
428 : integer :: spacecomm
429 : integer :: my_comm_atom, my_natom
430 : integer :: me
431 : integer :: natom
432 : ! integer :: ncid
433 : integer :: timrev
434 : character(len=fnlen) :: fname
435 : character(len=24) :: step_nb
436 : logical :: paral_atom
437 : logical :: remove_inv
438 : logical :: my_atmtab_allocated
439 69732 : type(crystal_t) :: crystal
440 1341 : type(ebands_t) :: ebands
441 : !arrays
442 1341 : integer, pointer :: my_atmtab(:)
443 1341 : real(dp), allocatable :: doccde(:)
444 :
445 : ! *************************************************************************
446 :
447 1341 : if (dtset%prtden /= 0) then
448 0 : spacecomm = mpi_enreg%comm_cell
449 0 : me = xmpi_comm_rank(spacecomm)
450 :
451 0 : natom = dtset%natom
452 0 : my_natom = mpi_enreg%my_natom
453 0 : paral_atom=(my_natom/=natom)
454 0 : my_comm_atom = mpi_enreg%comm_atom
455 0 : nullify(my_atmtab)
456 0 : if (paral_atom) then
457 0 : call get_my_atmtab(mpi_enreg%comm_atom,my_atmtab,my_atmtab_allocated,paral_atom,natom,my_natom_ref=my_natom)
458 : else
459 0 : ABI_MALLOC(my_atmtab, (natom))
460 0 : my_atmtab = (/ (iatom, iatom=1, natom) /)
461 0 : my_atmtab_allocated = .true.
462 : end if
463 :
464 : !FB: Maybe this should be moved out of that subroutine if needed in other outputs than densities
465 0 : remove_inv=.false.
466 0 : timrev = 2; if (any(dtset%kptopt == [3, 4])) timrev= 1
467 : call crystal%init(dtset%amu_orig(:,1),dtset%spgroup,natom,dtset%npsp,psps%ntypat, &
468 : dtset%nsym,tdks%rprimd,dtset%typat,tdks%xred,dtset%ziontypat,dtset%znucl,timrev,&
469 : dtset%nspden==2.and.dtset%nsppol==1,remove_inv,tdks%hdr%title,&
470 0 : dtset%symrel,dtset%tnons,dtset%symafm)
471 : !Electron band energies.
472 0 : bantot= dtset%mband*dtset%nkpt*dtset%nsppol
473 0 : ABI_CALLOC(doccde, (bantot))
474 : call ebands%init(bantot, dtset%nelect,dtset%ne_qFD,dtset%nh_qFD,dtset%ivalence, &
475 : doccde,tdks%eigen,dtset%istwfk,dtset%kptns,dtset%nband,dtset%nkpt,tdks%npwarr,dtset%nsppol, &
476 : dtset%nspinor,dtset%tphysel,dtset%tsmear,dtset%occopt,tdks%occ0,dtset%wtk,&
477 : dtset%cellcharge(1),dtset%kptopt,dtset%kptrlatt_orig,dtset%nshiftk_orig,dtset%shiftk_orig, &
478 0 : dtset%kptrlatt,dtset%nshiftk,dtset%shiftk)
479 0 : ABI_FREE(doccde)
480 :
481 0 : write(step_nb,*) istep
482 :
483 : !** Outputs the density
484 : !Warnings :
485 : !- core charge is excluded from the charge density;
486 : !- the potential is the INPUT vtrial.
487 0 : if (iwrite_fftdatar(mpi_enreg)) then
488 0 : fname = trim(dtfil%filnam_ds(4))//'_'//trim(adjustl(step_nb))//'_DEN'
489 : call fftdatar_write("density",fname,dtset%iomode,tdks%hdr,crystal,tdks%pawfgr%ngfft, &
490 0 : & cplex1,tdks%pawfgr%nfft,dtset%nspden,tdks%rhor,mpi_enreg,ebands=ebands)
491 : end if
492 :
493 0 : call crystal%free()
494 0 : call ebands%free()
495 0 : ABI_FREE(my_atmtab)
496 : end if
497 :
498 1341 : end subroutine prt_den
499 : !!***
500 :
501 : !!****f* m_rttddft_output/prt_dos
502 : !!
503 : !! NAME
504 : !! prt_dos
505 : !!
506 : !! FUNCTION
507 : !! Computes and outputs the electronic DOS
508 : !!
509 : !! INPUTS
510 : !! dtfil <type datafiles_type> = infos about file names, file unit numbers
511 : !! dtset <type(dataset_type)> = all input variables for this dataset
512 : !! istep <integer> = step number
513 : !! mpi_enreg <MPI_type> = MPI-parallelisation information
514 : !! psps <type(pseudopotential_type)> = variables related to pseudopotentials
515 : !! tdks <type(tdks_type)> = the tdks object to initialize
516 : !!
517 : !! OUTPUT
518 : !!
519 : !! SOURCE
520 1341 : subroutine prt_dos(dtfil, dtset, istep, mpi_enreg, psps, tdks)
521 :
522 : !Arguments ------------------------------------
523 : !scalars
524 : integer, intent(in) :: istep
525 : type(datafiles_type), intent(inout) :: dtfil
526 : type(dataset_type), intent(inout) :: dtset
527 : type(MPI_type), intent(inout) :: mpi_enreg
528 : type(pseudopotential_type), intent(inout) :: psps
529 : type(tdks_type), intent(inout) :: tdks
530 : !arrays
531 :
532 : !Local variables-------------------------------
533 : !scalars
534 : integer,parameter :: master=0
535 : integer :: bantot
536 : integer :: collect
537 : integer :: iatom
538 : integer :: spacecomm
539 : integer :: my_comm_atom, my_natom
540 : integer :: me
541 : integer :: natom
542 : ! integer :: ncid
543 : integer :: timrev
544 : character(len=fnlen) :: fname
545 : character(len=24) :: step_nb
546 : logical :: paral_atom
547 : logical :: remove_inv
548 : logical :: my_atmtab_allocated
549 75096 : type(crystal_t) :: crystal
550 1341 : type(epjdos_t) :: dos
551 1341 : type(ebands_t) :: ebands
552 : !arrays
553 1341 : integer, pointer :: my_atmtab(:)
554 1341 : real(dp), allocatable :: doccde(:)
555 :
556 : ! *************************************************************************
557 :
558 1341 : spacecomm = mpi_enreg%comm_cell
559 2682 : me = xmpi_comm_rank(spacecomm)
560 :
561 : !FB: @MT - Is this needed?
562 1341 : natom = dtset%natom
563 1341 : my_natom = mpi_enreg%my_natom
564 1341 : paral_atom=(my_natom/=natom)
565 1341 : my_comm_atom = mpi_enreg%comm_atom
566 1341 : nullify(my_atmtab)
567 1341 : if (paral_atom) then
568 0 : call get_my_atmtab(mpi_enreg%comm_atom,my_atmtab,my_atmtab_allocated,paral_atom,natom,my_natom_ref=my_natom)
569 : else
570 4023 : ABI_MALLOC(my_atmtab, (natom))
571 9006 : my_atmtab = (/ (iatom, iatom=1, natom) /)
572 1341 : my_atmtab_allocated = .true.
573 : end if
574 :
575 1341 : remove_inv=.false.
576 1341 : timrev = 2; if (any(dtset%kptopt == [3, 4])) timrev= 1
577 : call crystal%init(dtset%amu_orig(:,1),dtset%spgroup,natom,dtset%npsp,psps%ntypat, &
578 : dtset%nsym,tdks%rprimd,dtset%typat,tdks%xred,dtset%ziontypat,dtset%znucl,timrev,&
579 : dtset%nspden==2.and.dtset%nsppol==1,remove_inv,tdks%hdr%title,&
580 2682 : dtset%symrel,dtset%tnons,dtset%symafm)
581 : !Electron band energies.
582 1341 : bantot= dtset%mband*dtset%nkpt*dtset%nsppol
583 311767 : ABI_CALLOC(doccde, (bantot))
584 : call ebands%init(bantot,dtset%nelect,dtset%ne_qFD,dtset%nh_qFD,dtset%ivalence, &
585 : doccde,tdks%eigen,dtset%istwfk,dtset%kptns,dtset%nband,dtset%nkpt,tdks%npwarr,dtset%nsppol, &
586 : dtset%nspinor,dtset%tphysel,dtset%tsmear,dtset%occopt,tdks%occ0,dtset%wtk,&
587 : dtset%cellcharge(1),dtset%kptopt,dtset%kptrlatt_orig,dtset%nshiftk_orig,dtset%shiftk_orig, &
588 1341 : dtset%kptrlatt,dtset%nshiftk,dtset%shiftk)
589 1341 : ABI_FREE(doccde)
590 :
591 1341 : write(step_nb,*) istep
592 :
593 : !** Generate DOS using the tetrahedron method or using Gaussians
594 1341 : if (dtset%prtdos>=2.or.dtset%pawfatbnd>0) then
595 0 : call dos%init(dtset, psps, tdks%pawtab)
596 :
597 0 : if (dos%partial_dos_flag>=1 .or. dos%fatbands_flag==1)then
598 : ! Generate fractions for partial DOSs if needed partial_dos 1,2,3,4 give different decompositions
599 0 : collect = 1 !; if (psps%usepaw==1 .and. dos%partial_dos_flag /= 2) collect = 0
600 0 : if ((psps%usepaw==0.or.dtset%pawprtdos/=2) .and. dos%partial_dos_flag>=1) then
601 0 : call dos%partial_dos_fractions(crystal,dtset,tdks%eigen,tdks%occ0,tdks%npwarr,tdks%kg,tdks%cg,tdks%mcg,collect,mpi_enreg)
602 : end if
603 :
604 0 : if (psps%usepaw==1 .and. dos%partial_dos_flag /= 2) then
605 : ! TODO: update partial_dos_fractions_paw for extra atoms - no PAW contribution normally, but check bounds and so on.
606 0 : call dos%partial_dos_fractions_paw(tdks%cprj,tdks%dimcprj,dtset,tdks%mcprj,dtset%mkmem,mpi_enreg,tdks%pawrad,tdks%pawtab)
607 : end if
608 : else
609 0 : dos%fractions(:,:,:,1)=one
610 : end if
611 :
612 : !Here, print out fatbands for the k-points given in file appended _FATBANDS
613 0 : if (me == master .and. dtset%pawfatbnd>0 .and. dos%fatbands_flag==1) then
614 0 : fname = trim(dtfil%filnam_ds(4))//'_'//trim(adjustl(step_nb))//'_FATBANDS'
615 0 : call dos%prtfatbands(dtset,ebands,fname,dtset%pawfatbnd,tdks%pawtab)
616 : end if
617 :
618 : !Here, computation and output of DOS and partial DOS _DOS
619 0 : if (dos%fatbands_flag == 0 .and. dos%prtdos /= 4) then
620 0 : fname = trim(dtfil%filnam_ds(4))//'_'//trim(adjustl(step_nb))//'_DOS'
621 0 : call dos%calcnwrite(dtset,crystal,ebands,fname,spacecomm)
622 : end if
623 : end if
624 :
625 1341 : call dos%free()
626 1341 : call crystal%free()
627 1341 : call ebands%free()
628 1341 : ABI_FREE(my_atmtab)
629 :
630 2682 : end subroutine prt_dos
631 : !!***
632 :
633 : !!****f* m_rttddft_output/prt_wfk
634 : !!
635 : !! NAME
636 : !! prt_wfk
637 : !!
638 : !! FUNCTION
639 : !! Outputs wavefunctions in WFK file
640 : !!
641 : !! INPUTS
642 : !! dtfil <type datafiles_type> = infos about file names, file unit numbers
643 : !! dtset <type(dataset_type)> = all input variables for this dataset
644 : !! istep <integer> = step number
645 : !! mpi_enreg <MPI_type> = MPI-parallelisation information
646 : !! psps <type(pseudopotential_type)> = variables related to pseudopotentials
647 : !! tdks <type(tdks_type)> = the tdks object to initialize
648 : !! force_write <logical> = force the writing of WFK (useful for last step) - optional
649 : !!
650 : !! OUTPUT
651 : !!
652 : !! SOURCE
653 46 : subroutine prt_wfk(dtfil, dtset, istep, mpi_enreg, psps, tdks, force_write)
654 :
655 : !Arguments ------------------------------------
656 : !scalars
657 : integer, intent(in) :: istep
658 : type(datafiles_type), intent(inout) :: dtfil
659 : type(dataset_type), intent(inout) :: dtset
660 : type(MPI_type), intent(inout) :: mpi_enreg
661 : type(pseudopotential_type), intent(inout) :: psps
662 : type(tdks_type), intent(inout) :: tdks
663 : logical, optional, intent(in) :: force_write
664 : !arrays
665 :
666 : !Local variables-------------------------------
667 : !scalars
668 : integer,parameter :: response=0
669 : character(len=fnlen) :: fname
670 : character(len=24) :: step_nb
671 : logical :: lforce_write = .FALSE.
672 : !arrays
673 :
674 : ! *************************************************************************
675 :
676 46 : write(step_nb,*) istep
677 46 : fname = trim(dtfil%filnam_ds(4))//'_'//trim(adjustl(step_nb))//'_WFK'
678 :
679 46 : if (present(force_write)) then
680 46 : if (force_write) lforce_write = .TRUE.
681 : end if
682 :
683 : !Use initial eigenvalues to ensure that we get the same occupation upon restart
684 46 : if (lforce_write) then
685 : call outwf(tdks%cg,dtset,psps,tdks%eigen0,fname,tdks%hdr,tdks%kg,dtset%kptns, &
686 : & dtset%mband,tdks%mcg,dtset%mkmem,mpi_enreg,dtset%mpw,dtset%natom, &
687 : & dtset%nband,dtset%nkpt,tdks%npwarr,dtset%nsppol,tdks%occ0,response, &
688 46 : & dtfil%unwff2,tdks%wvl%wfs,tdks%wvl%descr, force_write=.TRUE.)
689 : else
690 : call outwf(tdks%cg,dtset,psps,tdks%eigen0,fname,tdks%hdr,tdks%kg,dtset%kptns, &
691 : & dtset%mband,tdks%mcg,dtset%mkmem,mpi_enreg,dtset%mpw,dtset%natom, &
692 : & dtset%nband,dtset%nkpt,tdks%npwarr,dtset%nsppol,tdks%occ0,response, &
693 0 : & dtfil%unwff2,tdks%wvl%wfs,tdks%wvl%descr)
694 : end if
695 :
696 46 : end subroutine prt_wfk
697 : !!***
698 :
699 : !!****f* m_rttddft_output/prt_restart
700 : !!
701 : !! NAME
702 : !! prt_restart
703 : !!
704 : !! FUNCTION
705 : !! Print restart file
706 : !!
707 : !! INPUTS
708 : !! dtfil <type datafiles_type> = infos about file names, file unit numbers
709 : !! istep <integer> = step number
710 : !! mpi_enreg <MPI_type> = MPI-parallelisation information
711 : !! tdks <type(tdks_type)> = the tdks object to initialize
712 : !!
713 : !! OUTPUT
714 : !!
715 : !! SOURCE
716 46 : subroutine prt_restart(dtfil, dtset, istep, mpi_enreg, tdks)
717 :
718 : !Arguments ------------------------------------
719 : !scalars
720 : integer, intent(in) :: istep
721 : type(datafiles_type), intent(inout) :: dtfil
722 : type(dataset_type), intent(inout) :: dtset
723 : type(MPI_type), intent(inout) :: mpi_enreg
724 : type(tdks_type), intent(inout) :: tdks
725 : !arrays
726 :
727 : !Local variables-------------------------------
728 : !scalars
729 : character(len=500) :: msg
730 : character(len=fnlen) :: fname
731 : character(len=24 ) :: step_nb
732 : !arrays
733 :
734 : ! *************************************************************************
735 :
736 46 : write(step_nb,*) istep
737 46 : if (mpi_enreg%me == 0) rewind(tdks%tdrestart_unit)
738 46 : write(msg,'(a)') step_nb
739 46 : call wrtout(tdks%tdrestart_unit,msg)
740 46 : write(msg,'(a)') trim(tdks%fname_wfk0)
741 46 : call wrtout(tdks%tdrestart_unit,msg)
742 46 : fname = trim(dtfil%filnam_ds(4))//'_'//trim(adjustl(step_nb))//'_WFK'
743 46 : write(msg,'(a)') trim(fname)
744 46 : call wrtout(tdks%tdrestart_unit,msg)
745 46 : write(msg,'(a)') trim(tdks%fname_tdener)
746 46 : call wrtout(tdks%tdrestart_unit,msg)
747 46 : if (dtset%td_ef_type /= 0) then
748 21 : write(msg,'(a)') trim(tdks%fname_tdef)
749 21 : call wrtout(tdks%tdrestart_unit,msg)
750 : end if
751 46 : if (dtset%prtcurrent /= 0) then
752 21 : write(msg,'(a)') trim(tdks%fname_current)
753 21 : call wrtout(tdks%tdrestart_unit,msg)
754 : end if
755 46 : if (dtset%td_ef_induced_vecpot /= 0) then
756 2 : write(msg,*) tdks%tdef%vecpot_ind(:,2)
757 2 : call wrtout(tdks%tdrestart_unit,msg)
758 2 : write(msg,*) tdks%tdef%vecpot_ind(:,1)
759 2 : call wrtout(tdks%tdrestart_unit,msg)
760 : end if
761 :
762 46 : end subroutine prt_restart
763 : !!***
764 :
765 : end module m_rttddft_output
766 : !!***
|
