Octopus
output.F90
Go to the documentation of this file.
1!! Copyright (C) 2002-2006 M. Marques, A. Castro, A. Rubio, G. Bertsch
2!!
3!! This program is free software; you can redistribute it and/or modify
4!! it under the terms of the GNU General Public License as published by
5!! the Free Software Foundation; either version 2, or (at your option)
6!! any later version.
7!!
8!! This program is distributed in the hope that it will be useful,
9!! but WITHOUT ANY WARRANTY; without even the implied warranty of
10!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
11!! GNU General Public License for more details.
12!!
13!! You should have received a copy of the GNU General Public License
14!! along with this program; if not, write to the Free Software
15!! Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
16!! 02110-1301, USA.
17!!
18
19#include "global.h"
20
22module output_oct_m
23 use accel_oct_m
24 use basins_oct_m
25 use box_oct_m
26 use comm_oct_m
27 use cube_oct_m
30 use debug_oct_m
33 use dos_oct_m
35 use elf_oct_m
36#if defined(HAVE_ETSF_IO)
37 use etsf_io
38 use etsf_io_tools
39#endif
42 use fft_oct_m
45 use global_oct_m
46 use grid_oct_m
49 use io_oct_m
51 use ions_oct_m
52 use, intrinsic :: iso_fortran_env
55 use lasers_oct_m
56 use lda_u_oct_m
59 use loct_oct_m
60 use math_oct_m
62 use mesh_oct_m
65 use mpi_oct_m
71 use parser_oct_m
74 use smear_oct_m
75 use space_oct_m
80 use stress_oct_m
81 use string_oct_m
85 use unit_oct_m
87 use utils_oct_m
89 use v_ks_oct_m
90 use vtk_oct_m
91 use xc_oct_m
92 use xc_oep_oct_m
94 use xc_f03_lib_m
95
96 implicit none
97
98 private
99 public :: &
100 output_bgw_t, &
101 output_init, &
104 output_all, &
106 doutput_lr, &
107 zoutput_lr, &
111
112contains
113
114 subroutine output_init(outp, namespace, space, st, gr, nst, ks)
115 type(output_t), intent(out) :: outp
116 type(namespace_t), intent(in) :: namespace
117 class(space_t), intent(in) :: space
118 type(states_elec_t), intent(in) :: st
119 type(grid_t), intent(in) :: gr
120 integer, intent(in) :: nst
121 type(v_ks_t), intent(inout) :: ks
122
123 type(block_t) :: blk
124 real(real64) :: norm
125 character(len=80) :: nst_string, default
126
127 push_sub(output_init)
128 outp%what = .false.
129
130 call io_function_read_what_how_when(namespace, space, outp%what, outp%how, outp%output_interval)
131
132 if (outp%what(option__output__wfs_fourier)) then
133 if (accel_is_enabled()) then
134 message(1) = "Wave functions in Fourier space not supported on GPUs."
135 call messages_fatal(1, namespace=namespace)
136 end if
137 call messages_experimental("Wave-functions in Fourier space", namespace=namespace)
138 end if
139
140 ! cannot calculate the ELF in 1D
141 if (outp%what(option__output__elf) .or. outp%what(option__output__elf_basins)) then
142 if (space%dim /= 2 .and. space%dim /= 3) then
143 outp%what(option__output__elf) = .false.
144 outp%what(option__output__elf_basins) = .false.
145 write(message(1), '(a)') 'Cannot calculate ELF except in 2D and 3D.'
146 call messages_warning(1, namespace=namespace)
147 end if
148 end if
149
150
151 if (outp%what(option__output__mmb_wfs)) then
152 call messages_experimental("Model many-body wfs", namespace=namespace)
153 end if
154
155 if (outp%what(option__output__xc_torque)) then
156 if (st%d%ispin /= spinors) then
157 write(message(1), '(a)') 'The output xc_torque can only be computed for spinors.'
158 call messages_fatal(1, namespace=namespace)
159 end if
160 if (space%dim /= 3) then
161 write(message(1), '(a)') 'The output xc_torque can only be computed in the 3D case.'
162 call messages_fatal(1, namespace=namespace)
163 end if
164 end if
165 if (outp%what(option__output__mmb_den)) then
166 call messages_experimental("Model many-body density matrix", namespace=namespace)
167 ! NOTES:
168 ! could be made into block to be able to specify which dimensions to trace
169 ! in principle all combinations are interesting, but this means we need to
170 ! be able to output density matrices for multiple particles or multiple
171 ! dimensions. The current 1D 1-particle case is simple.
172 end if
173
174 if (outp%what(option__output__energy_density)) call messages_experimental("'Output = energy_density'", namespace=namespace)
175 if (outp%what(option__output__heat_current)) call messages_experimental("'Output = heat_current'", namespace=namespace)
176
177 if (outp%what(option__output__wfs) .or. outp%what(option__output__wfs_sqmod)) then
178
179 !%Variable OutputWfsNumber
180 !%Type string
181 !%Default all states
182 !%Section Output
183 !%Description
184 !% Which wavefunctions to print, in list form: <i>i.e.</i>, "1-5" to print the first
185 !% five states, "2,3" to print the second and the third state, etc.
186 !% If more states are specified than available, extra ones will be ignored.
187 !%End
188
189 write(nst_string,'(i6)') nst
190 write(default,'(a,a)') "1-", trim(adjustl(nst_string))
191 call parse_variable(namespace, 'OutputWfsNumber', default, outp%wfs_list)
192 end if
193
194 if (parse_block(namespace, 'CurrentThroughPlane', blk) == 0) then
195 if (.not. outp%what(option__output__j_flow)) then
196 outp%what(option__output__j_flow) = .true.
197 call parse_variable(namespace, 'OutputInterval', 50, outp%output_interval(option__output__j_flow))
198 end if
199
200 !%Variable CurrentThroughPlane
201 !%Type block
202 !%Section Output
203 !%Description
204 !% The code can calculate current
205 !% traversing a user-defined portion of a plane, as specified by this block.
206 !% A small plain-text file <tt>current-flow</tt> will be written containing this information.
207 !% Only available for 1D, 2D, or 3D.
208 !% In the format below, <tt>origin</tt> is a point in the plane.
209 !% <tt>u</tt> and <tt>v</tt> are the (dimensionless) vectors defining the plane;
210 !% they will be normalized. <tt>spacing</tt> is the fineness of the mesh
211 !% on the plane. Integers <tt>nu</tt> and <tt>mu</tt> are the length and
212 !% width of the portion of the plane, in units of <tt>spacing</tt>.
213 !% Thus, the grid points included in the plane are
214 !% <tt>x_ij = origin + i*spacing*u + j*spacing*v</tt>,
215 !% for <tt>nu <= i <= mu</tt> and <tt>nv <= j <= mv</tt>.
216 !% Analogously, in the 2D case, the current flow is calculated through a line;
217 !% in the 1D case, the current flow is calculated through a point. Note that the spacing
218 !% can differ from the one used in the main calculation; an interpolation will be performed.
219 !%
220 !% Example (3D):
221 !%
222 !% <tt>%CurrentThroughPlane
223 !% <br>&nbsp;&nbsp; 0.0 | 0.0 | 0.0 # origin
224 !% <br>&nbsp;&nbsp; 0.0 | 1.0 | 0.0 # u
225 !% <br>&nbsp;&nbsp; 0.0 | 0.0 | 1.0 # v
226 !% <br>&nbsp;&nbsp; 0.2 # spacing
227 !% <br>&nbsp;&nbsp; 0 | 50 # nu | mu
228 !% <br>&nbsp;&nbsp; -50 | 50 # nv | mv
229 !% <br>%</tt>
230 !%
231 !% Example (2D):
232 !%
233 !% <tt>%CurrentThroughPlane
234 !% <br>&nbsp;&nbsp; 0.0 | 0.0 # origin
235 !% <br>&nbsp;&nbsp; 1.0 | 0.0 # u
236 !% <br>&nbsp;&nbsp; 0.2 # spacing
237 !% <br>&nbsp;&nbsp; 0 | 50 # nu | mu
238 !% <br>%</tt>
239 !%
240 !% Example (1D):
241 !%
242 !% <tt>%CurrentThroughPlane
243 !% <br>&nbsp;&nbsp; 0.0 # origin
244 !% <br>%</tt>
245 !%
246 !%End
247
248 select case (space%dim)
249 case (3)
250
251 call parse_block_float(blk, 0, 0, outp%plane%origin(1), units_inp%length)
252 call parse_block_float(blk, 0, 1, outp%plane%origin(2), units_inp%length)
253 call parse_block_float(blk, 0, 2, outp%plane%origin(3), units_inp%length)
254 call parse_block_float(blk, 1, 0, outp%plane%u(1))
255 call parse_block_float(blk, 1, 1, outp%plane%u(2))
256 call parse_block_float(blk, 1, 2, outp%plane%u(3))
257 call parse_block_float(blk, 2, 0, outp%plane%v(1))
258 call parse_block_float(blk, 2, 1, outp%plane%v(2))
259 call parse_block_float(blk, 2, 2, outp%plane%v(3))
260 call parse_block_float(blk, 3, 0, outp%plane%spacing, units_inp%length)
261 call parse_block_integer(blk, 4, 0, outp%plane%nu)
262 call parse_block_integer(blk, 4, 1, outp%plane%mu)
263 call parse_block_integer(blk, 5, 0, outp%plane%nv)
264 call parse_block_integer(blk, 5, 1, outp%plane%mv)
265
266 norm = norm2(outp%plane%u(1:3))
267 if (norm < m_epsilon) then
268 write(message(1), '(a)') 'u-vector for CurrentThroughPlane cannot have norm zero.'
269 call messages_fatal(1, namespace=namespace)
270 end if
271 outp%plane%u(1:3) = outp%plane%u(1:3) / norm
272
273 norm = norm2(outp%plane%v(1:3))
274 if (norm < m_epsilon) then
275 write(message(1), '(a)') 'v-vector for CurrentThroughPlane cannot have norm zero.'
276 call messages_fatal(1, namespace=namespace)
277 end if
278 outp%plane%v(1:3) = outp%plane%v(1:3) / norm
279
280 outp%plane%n(1) = outp%plane%u(2)*outp%plane%v(3) - outp%plane%u(3)*outp%plane%v(2)
281 outp%plane%n(2) = outp%plane%u(3)*outp%plane%v(1) - outp%plane%u(1)*outp%plane%v(3)
282 outp%plane%n(3) = outp%plane%u(1)*outp%plane%v(2) - outp%plane%u(2)*outp%plane%v(1)
283
284 case (2)
285
286 call parse_block_float(blk, 0, 0, outp%line%origin(1), units_inp%length)
287 call parse_block_float(blk, 0, 1, outp%line%origin(2), units_inp%length)
288 call parse_block_float(blk, 1, 0, outp%line%u(1))
289 call parse_block_float(blk, 1, 1, outp%line%u(2))
290 call parse_block_float(blk, 2, 0, outp%line%spacing, units_inp%length)
291 call parse_block_integer(blk, 3, 0, outp%line%nu)
292 call parse_block_integer(blk, 3, 1, outp%line%mu)
293
294 norm = norm2(outp%line%u(1:2))
295 if (norm < m_epsilon) then
296 write(message(1), '(a)') 'u-vector for CurrentThroughPlane cannot have norm zero.'
297 call messages_fatal(1, namespace=namespace)
298 end if
299 outp%line%u(1:2) = outp%line%u(1:2) / norm
300
301 outp%line%n(1) = -outp%line%u(2)
302 outp%line%n(2) = outp%line%u(1)
303
304 case (1)
305
306 call parse_block_float(blk, 0, 0, outp%line%origin(1), units_inp%length)
307
308 case default
309
310 call messages_not_implemented("CurrentThroughPlane for 4D or higher", namespace=namespace)
311
312 end select
313 call parse_block_end(blk)
314 end if
315
316 if (outp%what(option__output__matrix_elements)) then
317 call output_me_init(outp%me, namespace, space, st, gr, nst)
318 else
319 outp%me%what = .false.
320 end if
321
322 if (outp%what(option__output__berkeleygw)) then
323 if (accel_is_enabled()) then
324 message(1) = "BerkeleyGW is not compatible with GPUs."
325 call messages_fatal(1, namespace=namespace)
326 end if
327 call output_berkeleygw_init(nst, namespace, outp%bgw, space%periodic_dim)
328 end if
329
330 ! required for output_hamiltonian()
331 if (outp%what(option__output__potential_gradient) .and. .not. outp%what(option__output__potential)) then
332 outp%what(option__output__potential) = .true.
333 outp%output_interval(option__output__potential) = outp%output_interval(option__output__potential_gradient)
334 end if
335
336
337 !%Variable OutputDuringSCF
338 !%Type logical
339 !%Default no
340 !%Section Output
341 !%Description
342 !% During <tt>gs</tt> and <tt>unocc</tt> runs, if this variable is set to yes,
343 !% output will be written after every <tt>OutputInterval</tt> iterations.
344 !%End
345 call parse_variable(namespace, 'OutputDuringSCF', .false., outp%duringscf)
346
347 !%Variable RestartWriteInterval
348 !%Type integer
349 !%Default 50
350 !%Section Execution::IO
351 !%Description
352 !% Restart data is written when the iteration number is a multiple
353 !% of the <tt>RestartWriteInterval</tt> variable. For
354 !% time-dependent runs this includes the update of the output
355 !% controlled by the <tt>TDOutput</tt> variable. (Other output is
356 !% controlled by <tt>OutputInterval</tt>.)
357 !%End
358 call parse_variable(namespace, 'RestartWriteInterval', 50, outp%restart_write_interval)
359 if (outp%restart_write_interval <= 0) then
360 message(1) = "RestartWriteInterval must be > 0."
361 call messages_fatal(1, namespace=namespace)
362 end if
363
364 !%Variable OutputIterDir
365 !%Default "output_iter"
366 !%Type string
367 !%Section Output
368 !%Description
369 !% The name of the directory where <tt>Octopus</tt> stores information
370 !% such as the density, forces, etc. requested by variable <tt>Output</tt>
371 !% in the format specified by <tt>OutputFormat</tt>.
372 !% This information is written while iterating <tt>CalculationMode = gs</tt>, <tt>unocc</tt>, or <tt>td</tt>,
373 !% according to <tt>OutputInterval</tt>, and has nothing to do with the restart information.
374 !%End
375 call parse_variable(namespace, 'OutputIterDir', "output_iter", outp%iter_dir)
376 if (any(outp%what) .and. any(outp%output_interval > 0)) then
377 call io_mkdir(outp%iter_dir, namespace)
378 end if
379 call add_last_slash(outp%iter_dir)
380
381 ! At this point, we don`t know whether the states will be real or complex.
382 ! We therefore pass .false. to states_are_real, and need to check for real states later.
383
384 if (output_needs_current(outp, .false.)) then
385 call v_ks_calculate_current(ks, .true.)
386 else
387 call v_ks_calculate_current(ks, .false.)
388 end if
389
390 if (outp%what(option__output__current_dia)) then
391 message(1) = "The diamagnetic current will be calculated only if CalculateDiamagneticCurrent = yes."
392 call messages_warning(1, namespace=namespace)
393 end if
394
395 pop_sub(output_init)
396 end subroutine output_init
397
398 ! ---------------------------------------------------------
399 subroutine output_all(outp, namespace, space, dir, gr, ions, iter, st, hm, ks)
400 type(output_t), intent(in) :: outp
401 type(namespace_t), intent(in) :: namespace
402 class(space_t), intent(in) :: space
403 character(len=*), intent(in) :: dir
404 type(grid_t), intent(in) :: gr
405 type(ions_t), intent(in) :: ions
406 integer, intent(in) :: iter
407 type(states_elec_t), intent(inout) :: st
408 type(hamiltonian_elec_t), intent(inout) :: hm
409 type(v_ks_t), intent(inout) :: ks
410
411 integer :: idir, ierr, iout, iunit
412 character(len=MAX_PATH_LEN) :: fname
413
414 push_sub(output_all)
415 call profiling_in("OUTPUT_ALL")
416
417 if (any(outp%what)) then
418 message(1) = "Info: Writing output to " // trim(dir)
419 call messages_info(1, namespace=namespace)
420 call io_mkdir(dir, namespace)
421 end if
422
423 if (outp%what_now(option__output__mesh_r, iter)) then
424 do idir = 1, space%dim
425 write(fname, '(a,a)') 'mesh_r-', index2axis(idir)
426 call dio_function_output(outp%how(option__output__mesh_r), dir, fname, namespace, space, &
427 gr, gr%x(:,idir), units_out%length, ierr, pos=ions%pos, atoms=ions%atom)
428 end do
429 end if
430
431 call output_states(outp, namespace, space, dir, st, gr, ions, hm, iter)
432 call output_hamiltonian(outp, namespace, space, dir, hm, st, gr%der, ions, gr, iter, st%st_kpt_mpi_grp)
433
434 ! We can only test against the theory level here, as it is not set in the call of output_init().
435 if (outp%what_now(option__output__el_pressure, iter)) then
436 if(ks%theory_level /= kohn_sham_dft) then
437 call messages_not_implemented("el_pressure for TheoryLevel different from kohn_sham", namespace=namespace)
438 end if
439 if (st%d%spin_channels > 1) then
440 call messages_not_implemented("el_pressure for spin-polarized or spinors", namespace=namespace)
441 end if
442 end if
443
444 call output_localization_funct(outp, namespace, space, dir, st, hm, gr, ions, iter)
445
446 if (outp%what_now(option__output__j_flow, iter)) then
447 call output_current_flow(outp, namespace, space, dir, gr, st, hm%kpoints)
448 end if
449
450 if (outp%what_now(option__output__geometry, iter)) then
451 if (bitand(outp%how(option__output__geometry), option__outputformat__xcrysden) /= 0) then
452 call write_xsf_geometry_file(dir, "geometry", space, ions%latt, ions%pos, ions%atom, gr, namespace)
453 end if
454 if (bitand(outp%how(option__output__geometry), option__outputformat__xyz) /= 0) then
455 call ions%write_xyz(trim(dir)//'/geometry')
456 if (ions%space%is_periodic()) then
457 call ions%write_crystal(dir)
458 end if
459 end if
460 if (bitand(outp%how(option__output__geometry), option__outputformat__vtk) /= 0) then
461 call ions%write_vtk_geometry(trim(dir)//'/geometry')
462 end if
463 end if
464
465 if (outp%what_now(option__output__forces, iter)) then
466 if (bitand(outp%how(option__output__forces), option__outputformat__bild) /= 0) then
467 call ions%write_bild_forces_file(dir, "forces")
468 else
469 call write_xsf_geometry_file(dir, "forces", space, ions%latt, ions%pos, ions%atom, &
470 gr, namespace, total_forces = ions%tot_force)
471 end if
472 end if
473
474 if (outp%what_now(option__output__matrix_elements, iter)) then
475 call output_me(outp%me, namespace, space, dir, st, gr, ions, hm)
476 end if
477
478 do iout = lbound(outp%how, 1), ubound(outp%how, 1)
479 if (bitand(outp%how(iout), option__outputformat__etsf) /= 0) then
480 call output_etsf(outp, namespace, space, dir, st, gr, hm%kpoints, ions, iter)
481 exit
482 end if
483 end do
484
485 if (outp%what_now(option__output__berkeleygw, iter)) then
486 call output_berkeleygw(outp%bgw, namespace, space, dir, st, gr, ks, hm, ions)
487 end if
488
489 if (outp%what_now(option__output__energy_density, iter)) then
490 call output_energy_density(outp, namespace, space, dir, hm, ks, st, ions, gr)
491 end if
493 if (outp%what_now(option__output__stress, iter)) then
494 call io_mkdir(dir, namespace)
495 iunit = io_open(trim(dir)//'/stress', namespace, action='write')
496 call output_stress(iunit, space%periodic_dim, st%stress_tensors)
497 call io_close(iunit)
498 end if
499
500 if (hm%lda_u_level /= dft_u_none) then
501 if (outp%what_now(option__output__occ_matrices, iter))&
502 call lda_u_write_occupation_matrices(dir, hm%lda_u, st, namespace)
503
504 if (outp%what_now(option__output__effectiveu, iter))&
505 call lda_u_write_effectiveu(dir, hm%lda_u, namespace)
506
507 if (outp%what_now(option__output__magnetization, iter))&
508 call lda_u_write_magnetization(dir, hm%lda_u, ions, gr, st, namespace)
509
510 if (outp%what_now(option__output__local_orbitals, iter))&
511 call output_dftu_orbitals(outp, dir, namespace, space, hm%lda_u, st, gr, ions, hm%phase%is_allocated())
512
513 if (outp%what_now(option__output__kanamoriu, iter))&
514 call lda_u_write_kanamoriu(dir, st, hm%lda_u, namespace)
515
516 if (ks%oep_photon%level == oep_level_full) then
517 if (outp%what_now(option__output__photon_correlator, iter)) then
518 write(fname, '(a)') 'photon_correlator'
519 call dio_function_output(outp%how(option__output__photon_correlator), dir, trim(fname), namespace, space, &
520 gr, ks%oep_photon%pt%correlator(:,1), units_out%length, ierr, pos=ions%pos, atoms=ions%atom)
521 end if
522 end if
523 end if
524
525 ! We can only test against the theory level here, as it is not set in the call of output_init().
526 if (outp%what_now(option__output__xc_torque, iter)) then
527 if (ks%theory_level /= kohn_sham_dft .and. ks%theory_level /= generalized_kohn_sham_dft) then
528 write(message(1), '(a)') 'The output xc_torque can only be computed when there is a xc potential.'
529 call messages_fatal(1, namespace=namespace)
530 end if
531 end if
532
533 call output_xc_torque(outp, namespace, dir, gr, hm, st, ions, ions%space)
534
535 call profiling_out("OUTPUT_ALL")
536 pop_sub(output_all)
537 end subroutine output_all
538
539
540 ! ---------------------------------------------------------
541 subroutine output_localization_funct(outp, namespace, space, dir, st, hm, gr, ions, iter)
542 type(output_t), intent(in) :: outp
543 type(namespace_t), intent(in) :: namespace
544 class(space_t), intent(in) :: space
545 character(len=*), intent(in) :: dir
546 type(states_elec_t), intent(inout) :: st
547 type(hamiltonian_elec_t), intent(in) :: hm
548 type(grid_t), intent(in) :: gr
549 type(ions_t), intent(in) :: ions
550 integer, intent(in) :: iter
551
552 real(real64), allocatable :: f_loc(:,:)
553 character(len=MAX_PATH_LEN) :: fname
554 integer :: is, ierr, imax
555 type(mpi_grp_t) :: mpi_grp
556
558
559 mpi_grp = st%dom_st_kpt_mpi_grp
560
561 ! if SPIN_POLARIZED, the ELF contains one extra channel: the total ELF
562 imax = st%d%nspin
563 if (st%d%ispin == spin_polarized) imax = 3
564
565 safe_allocate(f_loc(1:gr%np, 1:imax))
566
567 ! First the ELF in real space
568 if (outp%what_now(option__output__elf, iter) .or. outp%what_now(option__output__elf_basins, iter)) then
569 assert(space%dim /= 1)
570
571 call elf_calc(space, st, gr, hm%kpoints, f_loc)
572
573 ! output ELF in real space
574 if (outp%what_now(option__output__elf, iter)) then
575 write(fname, '(a)') 'elf_rs'
576 call dio_function_output(outp%how(option__output__elf), dir, trim(fname), namespace, space, gr, &
577 f_loc(:,imax), unit_one, ierr, pos=ions%pos, atoms=ions%atom, grp = mpi_grp)
578 ! this quantity is dimensionless
579
580 if (st%d%ispin /= unpolarized) then
581 do is = 1, 2
582 write(fname, '(a,i1)') 'elf_rs-sp', is
583 call dio_function_output(outp%how(option__output__elf), dir, trim(fname), namespace, space, gr, &
584 f_loc(:, is), unit_one, ierr, pos=ions%pos, atoms=ions%atom, grp = mpi_grp)
585 ! this quantity is dimensionless
586 end do
587 end if
588 end if
589
590 if (outp%what_now(option__output__elf_basins, iter)) then
591 call out_basins(f_loc(:,1), "elf_rs_basins", outp%how(option__output__elf_basins))
592 end if
593 end if
594
595 ! Now Bader analysis
596 if (outp%what_now(option__output__bader, iter)) then
597 do is = 1, st%d%nspin
598 call dderivatives_lapl(gr%der, st%rho(:,is), f_loc(:,is))
599
600 fname = get_filename_with_spin('bader', st%d%nspin, is)
601
602 call dio_function_output(outp%how(option__output__bader), dir, trim(fname), namespace, space, gr, &
603 f_loc(:,is), units_out%length**(-2 - space%dim), ierr, &
604 pos=ions%pos, atoms=ions%atom, grp = mpi_grp)
605
606 fname = get_filename_with_spin('bader_basins', st%d%nspin, is)
607 call out_basins(f_loc(:,1), fname, outp%how(option__output__bader))
608 end do
609 end if
610
611 ! Now the pressure
612 if (outp%what_now(option__output__el_pressure, iter)) then
613 call calc_electronic_pressure(st, hm, gr, f_loc(:,1))
614 call dio_function_output(outp%how(option__output__el_pressure), dir, "el_pressure", namespace, space, gr, &
615 f_loc(:,1), unit_one, ierr, pos=ions%pos, atoms=ions%atom, grp = mpi_grp)
616 ! this quantity is dimensionless
617 end if
618
619 safe_deallocate_a(f_loc)
620
622
623 contains
624 ! ---------------------------------------------------------
625 subroutine out_basins(ff, filename, output_how)
626 real(real64), intent(in) :: ff(:)
627 character(len=*), intent(in) :: filename
628 integer(int64), intent(in) :: output_how
629
630 character(len=MAX_PATH_LEN) :: fname
631 type(basins_t) :: basins
632 integer :: iunit
633
635
636 call basins_init(basins, namespace, gr)
637 call basins_analyze(basins, namespace, gr, ff(:), st%rho, 0.01_real64)
638
639 call dio_function_output(output_how, dir, trim(filename), namespace, space, gr, &
640 real(basins%map, real64) , unit_one, ierr, pos=ions%pos, atoms=ions%atom, grp = mpi_grp)
641 ! this quantity is dimensionless
642
643 write(fname,'(4a)') trim(dir), '/', trim(filename), '.info'
644 iunit = io_open(trim(fname), namespace, action = 'write')
645 call basins_write(basins, gr, iunit)
646 call io_close(iunit)
647
648 call basins_end(basins)
649
651 end subroutine out_basins
652
653 end subroutine output_localization_funct
654
655
656 ! ---------------------------------------------------------
657 subroutine calc_electronic_pressure(st, hm, gr, pressure)
658 type(states_elec_t), intent(inout) :: st
659 type(hamiltonian_elec_t), intent(in) :: hm
660 type(grid_t), intent(in) :: gr
661 real(real64), intent(out) :: pressure(:)
662
663 real(real64), allocatable :: rho(:,:), lrho(:), tau(:,:)
664 real(real64) :: p_tf, dens
665 integer :: is, ii
666
668
669 safe_allocate( rho(1:gr%np_part, 1:st%d%nspin))
670 safe_allocate(lrho(1:gr%np))
671 safe_allocate( tau(1:gr%np, 1:st%d%nspin))
672
673 rho = m_zero
674 call density_calc(st, gr, rho)
675 call states_elec_calc_quantities(gr, st, hm%kpoints, .false., kinetic_energy_density = tau)
676
677 pressure = m_zero
678 do is = 1, st%d%spin_channels
679 lrho = m_zero
680 call dderivatives_lapl(gr%der, rho(:, is), lrho)
681
682 pressure(:) = pressure(:) + &
683 tau(:, is)/m_three - lrho(:)/m_four
684 end do
685
686 do ii = 1, gr%np
687 dens = sum(rho(ii,1:st%d%spin_channels))
688
689 p_tf = m_two/m_five*(m_three*m_pi**2)**(m_two/m_three)* &
690 dens**(m_five/m_three)
691
692 ! add XC pressure
693 ! FIXME: Not correct for spinors and spin-polarized here
694 pressure(ii) = pressure(ii) + (dens*hm%ks_pot%vxc(ii,1) - hm%energy%exchange - hm%energy%correlation)
695
696 pressure(ii) = pressure(ii)/p_tf
697 pressure(ii) = m_half*(m_one + pressure(ii)/sqrt(m_one + pressure(ii)**2))
698 end do
699
701 end subroutine calc_electronic_pressure
702
703
704 ! ---------------------------------------------------------
705 subroutine output_energy_density(outp, namespace, space, dir, hm, ks, st, ions, gr)
706 type(output_t), intent(in) :: outp
707 type(namespace_t), intent(in) :: namespace
708 class(space_t), intent(in) :: space
709 character(len=*), intent(in) :: dir
710 type(hamiltonian_elec_t), intent(in) :: hm
711 type(v_ks_t), intent(inout) :: ks
712 type(states_elec_t), intent(in) :: st
713 type(ions_t), intent(in) :: ions
714 type(grid_t), intent(in) :: gr
715
716 integer :: is, ierr, ip
717 character(len=MAX_PATH_LEN) :: fname
718 type(unit_t) :: fn_unit
719 real(real64), allocatable :: energy_density(:, :)
720 real(real64), allocatable :: ex_density(:)
721 real(real64), allocatable :: ec_density(:)
722
723 push_sub(output_energy_density)
724
725 fn_unit = units_out%energy*units_out%length**(-space%dim)
726 safe_allocate(energy_density(1:gr%np, 1:st%d%nspin))
727
728 ! the kinetic energy density
729 call states_elec_calc_quantities(gr, st, hm%kpoints, .true., kinetic_energy_density = energy_density)
730
731 ! the external potential energy density
732 do is = 1, st%d%nspin
733 do ip = 1, gr%np
734 energy_density(ip, is) = energy_density(ip, is) + st%rho(ip, is)*hm%ep%vpsl(ip)
735 end do
736 end do
737
738 ! the hartree energy density
739 do is = 1, st%d%nspin
740 do ip = 1, gr%np
741 energy_density(ip, is) = energy_density(ip, is) + m_half*st%rho(ip, is)*hm%ks_pot%vhartree(ip)
742 end do
743 end do
744
745 ! the XC energy density
746 safe_allocate(ex_density(1:gr%np))
747 safe_allocate(ec_density(1:gr%np))
748
749 call xc_get_vxc(gr, ks%xc, st, hm%kpoints, hm%psolver, namespace, space, st%rho, st%d%ispin, &
750 hm%ions%latt%rcell_volume, ex_density = ex_density, ec_density = ec_density)
751 do is = 1, st%d%nspin
752 do ip = 1, gr%np
753 energy_density(ip, is) = energy_density(ip, is) + ex_density(ip) + ec_density(ip)
754 end do
755 end do
756
757 safe_deallocate_a(ex_density)
758 safe_deallocate_a(ec_density)
759
760 do is = 1, st%d%spin_channels
761 fname = get_filename_with_spin('energy_density', st%d%nspin, is)
762 call dio_function_output(outp%how(option__output__energy_density), dir, trim(fname), namespace, space, gr, &
763 energy_density(:, is), unit_one, ierr, pos=ions%pos, atoms=ions%atom, grp = st%dom_st_kpt_mpi_grp)
764 end do
765 safe_deallocate_a(energy_density)
766
767 pop_sub(output_energy_density)
768 end subroutine output_energy_density
769
770 !--------------------------------------------------------------
771
772 logical function output_need_exchange(outp) result(need_exx)
773 type(output_t), intent(in) :: outp
774
775 need_exx =(outp%what(option__output__berkeleygw) &
776 .or. outp%me%what(option__outputmatrixelements__two_body) &
777 .or. outp%me%what(option__outputmatrixelements__two_body_exc_k))
778 end function output_need_exchange
779
780
781 ! ---------------------------------------------------------
782 subroutine output_dftu_orbitals(outp, dir, namespace, space, this, st, mesh, ions, has_phase)
783 type(output_t), intent(in) :: outp
784 character(len=*), intent(in) :: dir
785 type(namespace_t), intent(in) :: namespace
786 class(space_t), intent(in) :: space
787 type(lda_u_t), intent(in) :: this
788 type(states_elec_t), intent(in) :: st
789 class(mesh_t), intent(in) :: mesh
790 type(ions_t), intent(in) :: ions
791 logical, intent(in) :: has_phase
792
793 integer :: ios, im, ik, idim, ierr
794 complex(real64), allocatable :: tmp(:)
795 real(real64), allocatable :: dtmp(:)
796 type(orbitalset_t), pointer :: os
797 type(unit_t) :: fn_unit
798 character(len=MAX_PATH_LEN) :: fname
799
800 push_sub(output_dftu_orbitals)
801
802 fn_unit = sqrt(units_out%length**(-space%dim))
803
804 if (this%basis%use_submesh) then
805 if (states_are_real(st)) then
806 safe_allocate(dtmp(1:mesh%np))
807 else
808 safe_allocate(tmp(1:mesh%np))
809 end if
810 end if
811
812 do ios = 1, this%norbsets
813 os => this%orbsets(ios)
814 do ik = st%d%kpt%start, st%d%kpt%end
815 do im = 1, this%orbsets(ios)%norbs
816 do idim = 1, min(os%ndim, st%d%dim)
817 if (st%nik > 1) then
818 if (min(os%ndim, st%d%dim) > 1) then
819 write(fname, '(a,i1,a,i3.3,a,i8.8,a,i1)') 'orb', im, '-os', ios, '-k', ik, '-sp', idim
820 else
821 write(fname, '(a,i1,a,i3.3,a,i8.8)') 'orb', im, '-os', ios, '-k', ik
822 end if
823 else
824 if (min(os%ndim, st%d%dim) > 1) then
825 write(fname, '(a,i1,a,i3.3,a,i1)') 'orb', im, '-os', ios, '-sp', idim
826 else
827 write(fname, '(a,i1,a,i3.3)') 'orb', im, '-os', ios
828 end if
829 end if
830 if (has_phase) then
831 if (.not. this%basis%use_submesh) then
832 call zio_function_output(outp%how(option__output__local_orbitals), dir, fname, namespace, space, &
833 mesh, os%eorb_mesh(1:mesh%np,im,idim,ik), fn_unit, ierr, pos=ions%pos, atoms=ions%atom)
834 else
835 tmp = m_z0
836 call submesh_add_to_mesh(os%sphere, os%eorb_submesh(1:os%sphere%np,idim,im,ik), tmp)
837 call zio_function_output(outp%how(option__output__local_orbitals), dir, fname, namespace, space, &
838 mesh, tmp, fn_unit, ierr, pos=ions%pos, atoms=ions%atom)
839 end if
840 else
841 if (.not. this%basis%use_submesh) then
842 if (states_are_real(st)) then
843 call dio_function_output(outp%how(option__output__local_orbitals), dir, fname, namespace, space, mesh, &
844 os%dorb(1:mesh%np,idim,im), fn_unit, ierr, pos=ions%pos, atoms=ions%atom)
845 else
846 call zio_function_output(outp%how(option__output__local_orbitals), dir, fname, namespace, space, mesh, &
847 os%zorb(1:mesh%np,idim,im), fn_unit, ierr, pos=ions%pos, atoms=ions%atom)
848 end if
849 else
850 if (states_are_real(st)) then
851 dtmp = m_z0
852 call submesh_add_to_mesh(os%sphere, os%dorb(1:os%sphere%np,idim,im), dtmp)
853 call dio_function_output(outp%how(option__output__local_orbitals), dir, fname, namespace, space, &
854 mesh, dtmp, fn_unit, ierr, pos=ions%pos, atoms=ions%atom)
855 else
856 tmp = m_z0
857 call submesh_add_to_mesh(os%sphere, os%zorb(1:os%sphere%np,idim,im), tmp)
858 call zio_function_output(outp%how(option__output__local_orbitals), dir, fname, namespace, space, &
859 mesh, tmp, fn_unit, ierr, pos=ions%pos, atoms=ions%atom)
860 end if
861 end if
862 end if
863 end do
864 end do
865 end do
866 end do
867
868 safe_deallocate_a(tmp)
869 safe_deallocate_a(dtmp)
870
871 pop_sub(output_dftu_orbitals)
872 end subroutine output_dftu_orbitals
873
874 ! ---------------------------------------------------------
875 logical function output_needs_current(outp, states_are_real)
876 type(output_t), intent(in) :: outp
877 logical, intent(in) :: states_are_real
878
879 output_needs_current = .false.
880
881 if (outp%what(option__output__current) &
882 .or. outp%what(option__output__current_dia) &
883 .or. outp%what(option__output__heat_current) &
884 .or. outp%what(option__output__current_kpt)) then
885 if (.not. states_are_real) then
887 else
888 message(1) = 'No current density output for real states since it is identically zero.'
889 call messages_warning(1)
890 end if
891 end if
892
893
894 end function
895
896#include "output_etsf_inc.F90"
897
898#include "output_states_inc.F90"
899
900#include "output_h_inc.F90"
901
902#include "undef.F90"
903#include "complex.F90"
904#include "output_linear_response_inc.F90"
905
906#include "undef.F90"
907#include "real.F90"
908#include "output_linear_response_inc.F90"
909
910end module output_oct_m
911
912!! Local Variables:
913!! mode: f90
914!! coding: utf-8
915!! End:
pure logical function, public accel_is_enabled()
Definition: accel.F90:400
subroutine, public basins_write(this, mesh, iunit)
Definition: basins.F90:352
subroutine, public basins_init(this, namespace, mesh)
Definition: basins.F90:151
subroutine, public basins_analyze(this, namespace, mesh, f, rho, threshold)
Definition: basins.F90:188
subroutine, public basins_end(this)
Definition: basins.F90:170
This module implements a calculator for the density and defines related functions.
Definition: density.F90:120
subroutine, public density_calc(st, gr, density, istin)
Computes the density from the orbitals in st.
Definition: density.F90:609
This module calculates the derivatives (gradients, Laplacians, etc.) of a function.
subroutine, public dderivatives_lapl(der, ff, op_ff, ghost_update, set_bc, factor)
apply the Laplacian to a mesh function
Module that handles computing and output of various density of states.
Definition: dos.F90:116
integer, parameter, public unpolarized
Parameters...
integer, parameter, public spinors
integer, parameter, public spin_polarized
subroutine, public elf_calc(space, st, gr, kpoints, elf, de)
(time-dependent) electron localization function, (TD)ELF.
Definition: elf.F90:167
Fast Fourier Transform module. This module provides a single interface that works with different FFT ...
Definition: fft.F90:118
real(real64), parameter, public m_two
Definition: global.F90:190
real(real64), parameter, public m_zero
Definition: global.F90:188
real(real64), parameter, public m_four
Definition: global.F90:192
real(real64), parameter, public m_pi
some mathematical constants
Definition: global.F90:186
complex(real64), parameter, public m_z0
Definition: global.F90:198
real(real64), parameter, public m_epsilon
Definition: global.F90:204
real(real64), parameter, public m_half
Definition: global.F90:194
real(real64), parameter, public m_one
Definition: global.F90:189
real(real64), parameter, public m_three
Definition: global.F90:191
real(real64), parameter, public m_five
Definition: global.F90:193
This module implements the underlying real-space grid.
Definition: grid.F90:117
This module defines classes and functions for interaction partners.
subroutine, public zio_function_output(how, dir, fname, namespace, space, mesh, ff, unit, ierr, pos, atoms, grp, root)
Top-level IO routine for functions defined on the mesh.
subroutine, public io_function_read_what_how_when(namespace, space, what, how, output_interval, what_tag_in, how_tag_in, output_interval_tag_in, ignore_error)
subroutine, public dio_function_output(how, dir, fname, namespace, space, mesh, ff, unit, ierr, pos, atoms, grp, root)
Top-level IO routine for functions defined on the mesh.
subroutine, public write_xsf_geometry_file(dir, fname, space, latt, pos, atoms, mesh, namespace, total_forces)
Definition: io.F90:114
subroutine, public io_close(iunit, grp)
Definition: io.F90:418
subroutine, public io_mkdir(fname, namespace, parents)
Definition: io.F90:311
integer function, public io_open(file, namespace, action, status, form, position, die, recl, grp)
Definition: io.F90:352
A module to handle KS potential, without the external potential.
integer, parameter, public generalized_kohn_sham_dft
integer, parameter, public kohn_sham_dft
subroutine, public lda_u_write_occupation_matrices(dir, this, st, namespace)
Prints the occupation matrices at the end of the scf calculation.
Definition: lda_u_io.F90:157
subroutine, public lda_u_write_kanamoriu(dir, st, this, namespace)
Definition: lda_u_io.F90:352
subroutine, public lda_u_write_magnetization(dir, this, ions, mesh, st, namespace)
Definition: lda_u_io.F90:483
subroutine, public lda_u_write_effectiveu(dir, this, namespace)
Definition: lda_u_io.F90:265
integer, parameter, public dft_u_none
Definition: lda_u.F90:201
This module is intended to contain "only mathematical" functions and procedures.
Definition: math.F90:115
This module defines various routines, operating on mesh functions.
This module defines the meshes, which are used in Octopus.
Definition: mesh.F90:118
subroutine, public messages_not_implemented(feature, namespace)
Definition: messages.F90:1113
subroutine, public messages_warning(no_lines, all_nodes, namespace)
Definition: messages.F90:537
character(len=256), dimension(max_lines), public message
to be output by fatal, warning
Definition: messages.F90:160
subroutine, public messages_fatal(no_lines, only_root_writes, namespace)
Definition: messages.F90:414
subroutine, public messages_experimental(name, namespace)
Definition: messages.F90:1085
subroutine, public messages_info(no_lines, iunit, debug_only, stress, all_nodes, namespace)
Definition: messages.F90:616
subroutine, public output_berkeleygw(bgw, namespace, space, dir, st, gr, ks, hm, ions)
subroutine, public output_berkeleygw_init(nst, namespace, bgw, periodic_dim)
this module contains the low-level part of the output system
Definition: output_low.F90:115
character(len=max_path_len) function, public get_filename_with_spin(output, nspin, spin_index)
Returns the filame as output, or output-spX is spin polarized.
Definition: output_low.F90:232
subroutine, public output_me_init(this, namespace, space, st, gr, nst)
Definition: output_me.F90:155
subroutine, public output_me(this, namespace, space, dir, st, gr, ions, hm)
Definition: output_me.F90:278
this module contains the output system
Definition: output.F90:115
subroutine calc_electronic_pressure(st, hm, gr, pressure)
Definition: output.F90:751
subroutine, public output_states(outp, namespace, space, dir, st, gr, ions, hm, iter)
Definition: output.F90:1888
logical function, public output_needs_current(outp, states_are_real)
Definition: output.F90:969
subroutine, public output_hamiltonian(outp, namespace, space, dir, hm, st, der, ions, gr, iter, grp)
Definition: output.F90:2184
subroutine, public output_all(outp, namespace, space, dir, gr, ions, iter, st, hm, ks)
Definition: output.F90:493
logical function, public output_need_exchange(outp)
Definition: output.F90:866
subroutine, public output_init(outp, namespace, space, st, gr, nst, ks)
Definition: output.F90:208
subroutine output_xc_torque(outp, namespace, dir, mesh, hm, st, ions, space)
Definition: output.F90:2502
subroutine, public output_current_flow(outp, namespace, space, dir, gr, st, kpoints)
Definition: output.F90:2076
subroutine, public zoutput_lr(outp, namespace, space, dir, st, mesh, lr, idir, isigma, ions, pert_unit)
Definition: output.F90:2609
subroutine, public doutput_lr(outp, namespace, space, dir, st, mesh, lr, idir, isigma, ions, pert_unit)
Definition: output.F90:2878
subroutine, public output_scalar_pot(outp, namespace, space, dir, mesh, ions, ext_partners, time)
Definition: output.F90:2465
subroutine output_energy_density(outp, namespace, space, dir, hm, ks, st, ions, gr)
Definition: output.F90:799
subroutine output_dftu_orbitals(outp, dir, namespace, space, this, st, mesh, ions, has_phase)
Definition: output.F90:876
subroutine output_localization_funct(outp, namespace, space, dir, st, hm, gr, ions, iter)
Definition: output.F90:635
subroutine output_etsf(outp, namespace, space, dir, st, gr, kpoints, ions, iter)
To create an etsf file one has to do the following:
Definition: output.F90:1020
integer function, public parse_block(namespace, name, blk, check_varinfo_)
Definition: parser.F90:618
subroutine, public profiling_out(label)
Increment out counter and sum up difference between entry and exit time.
Definition: profiling.F90:623
subroutine, public profiling_in(label, exclude)
Increment in counter and save entry time.
Definition: profiling.F90:552
pure logical function, public states_are_real(st)
This module defines routines to write information about states.
subroutine, public states_elec_calc_quantities(gr, st, kpoints, nlcc, kinetic_energy_density, paramagnetic_current, density_gradient, density_laplacian, gi_kinetic_energy_density, st_end)
calculated selected quantities
This module implements the calculation of the stress tensor.
Definition: stress.F90:118
subroutine, public output_stress(iunit, space_dim, stress_tensors, all_terms)
Definition: stress.F90:1051
subroutine, public add_last_slash(str)
Adds a '/' in the end of the string, only if it missing. Useful for directories.
Definition: string.F90:161
brief This module defines the class unit_t which is used by the unit_systems_oct_m module.
Definition: unit.F90:132
This module defines the unit system, used for input and output.
type(unit_system_t), public units_out
type(unit_system_t), public units_inp
the units systems for reading and writing
type(unit_t), public unit_one
some special units required for particular quantities
This module is intended to contain simple general-purpose utility functions and procedures.
Definition: utils.F90:118
character pure function, public index2axis(idir)
Definition: utils.F90:202
subroutine, public v_ks_calculate_current(this, calc_cur)
Definition: v_ks.F90:1436
integer, parameter, public oep_level_full
Definition: xc_oep.F90:172
subroutine out_basins(ff, filename, output_how)
Definition: output.F90:719
Description of the grid, containing information on derivatives, stencil, and symmetries.
Definition: grid.F90:168
Class to describe DFT+U parameters.
Definition: lda_u.F90:214
Describes mesh distribution to nodes.
Definition: mesh.F90:186
This is defined even when running serial.
Definition: mpi.F90:142
Output information for BerkeleyGW.
Definition: output_low.F90:145
output handler class
Definition: output_low.F90:164
The states_elec_t class contains all electronic wave functions.
int true(void)