Octopus
td.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
21module td_oct_m
27 use debug_oct_m
33 use epot_oct_m
35 use forces_oct_m
37 use global_oct_m
38 use grid_oct_m
41 use io_oct_m
43 use ions_oct_m
44 use kick_oct_m
45 use, intrinsic :: iso_fortran_env
46 use lasers_oct_m
47 use lda_u_oct_m
50 use loct_oct_m
52 use mesh_oct_m
54 use mpi_oct_m
57 use output_oct_m
59 use parser_oct_m
60 use pes_oct_m
70 use scf_oct_m
72 use space_oct_m
76 use stress_oct_m
78 use types_oct_m
79 use unit_oct_m
81 use v_ks_oct_m
84 use xc_oct_m
85
86 implicit none
87
88 private
89 public :: &
90 td_t, &
91 td_run, &
92 td_init, &
94 td_end, &
95 td_end_run, &
98 td_dump, &
106
108 integer, parameter, public :: &
109 EHRENFEST = 1, &
110 bo = 2
111
112 type td_t
113 private
114 type(propagator_base_t), public :: tr
115 type(scf_t), public :: scf
116 type(ion_dynamics_t), public :: ions_dyn
117 real(real64), public :: dt
118 integer, public :: max_iter
119 integer, public :: iter
120 logical, public :: recalculate_gs
121
122 type(pes_t), public :: pesv
123
124 integer, public :: dynamics
125 integer, public :: energy_update_iter
126 real(real64) :: scissor
127
128 logical :: freeze_occ
129 logical :: freeze_u
130 integer :: freeze_orbitals
131
132 logical :: from_scratch = .false.
133
134 type(td_write_t), public :: write_handler
135 type(restart_t) :: restart_load
136 type(restart_t) :: restart_dump
137
138 type(dmp_t) :: dmp
139 end type td_t
140
141
142contains
143
144 subroutine td_init(td, namespace, space, gr, ions, st, ks, hm, ext_partners, outp)
145 type(td_t), intent(inout) :: td
146 type(namespace_t), intent(in) :: namespace
147 class(space_t), intent(in) :: space
148 type(grid_t), intent(in) :: gr
149 type(ions_t), intent(inout) :: ions
150 type(states_elec_t), intent(in) :: st
151 type(v_ks_t), intent(in) :: ks
152 type(hamiltonian_elec_t), intent(in) :: hm
153 type(partner_list_t), intent(in) :: ext_partners
154 type(output_t), intent(in) :: outp
155
156 integer :: default
157 real(real64) :: propagation_time
158 type(lasers_t), pointer :: lasers
159 logical :: symmetrize
160
161 push_sub(td_init)
162
163 if (hm%pcm%run_pcm) call messages_experimental("PCM for CalculationMode = td", namespace=namespace)
164
165 symmetrize = hm%kpoints%use_symmetries .or. st%symmetrize_density
166 call ion_dynamics_init(td%ions_dyn, namespace, ions, symmetrize, gr%symmetrizer)
167
168 if (td%ions_dyn%ions_move()) then
169 if (hm%kpoints%use_symmetries) then
170 message(1) = "KPoints symmetries cannot be used with moving ions."
171 message(2) = "Please set KPointsSymmetries = no."
172 call messages_fatal(2, namespace=namespace)
173 end if
174 if (st%symmetrize_density) then
175 message(1) = "Symmetrization of the density cannot be used with moving ions."
176 message(2) = "Please set SymmetrizeDensity = no."
177 call messages_fatal(2, namespace=namespace)
178 end if
179 end if
180
181 td%iter = 0
182
183 !%Variable TDTimeStep
184 !%Type float
185 !%Section Time-Dependent::Propagation
186 !%Description
187 !% The time-step for the time propagation. For most propagators you
188 !% want to use the largest value that is possible without the
189 !% evolution becoming unstable.
190 !%
191 !% While prior versions of Octopus used to have a default time step, however,
192 !% this now needs to be systematically defined.
193 !%End
194
195 call parse_variable(namespace, 'TDTimeStep', -m_one, td%dt, unit = units_inp%time)
196
197 if (td%dt <= m_zero) then
198 write(message(1),'(a)') 'A positive value for TDTimeStep must be defined in the input file.'
199 call messages_fatal(1, namespace=namespace)
200 end if
201
202 call messages_print_var_value('TDTimeStep', td%dt, unit = units_out%time, namespace=namespace)
204
205 if (parse_is_defined(namespace, 'TDMaxSteps') .and. parse_is_defined(namespace, 'TDPropagationTime')) then
206 call messages_write('You cannot set TDMaxSteps and TDPropagationTime at the same time')
207 call messages_fatal(namespace=namespace)
208 end if
210 !%Variable TDPropagationTime
211 !%Type float
212 !%Section Time-Dependent::Propagation
213 !%Description
214 !% The length of the time propagation. You cannot set this variable
215 !% at the same time as <tt>TDMaxSteps</tt>. By default this variable will
216 !% not be used.
217 !%
218 !% The units for this variable are <math>\hbar</math>/Hartree (or <math>\hbar</math>/eV if you
219 !% selected <tt>ev_angstrom</tt> as input units). The approximate conversions to
220 !% femtoseconds are 1 fs = 41.34 <math>\hbar</math>/Hartree = 1.52 <math>\hbar</math>/eV.
221 !%End
222 call parse_variable(namespace, 'TDPropagationTime', -1.0_real64, propagation_time, unit = units_inp%time)
224 call messages_obsolete_variable(namespace, 'TDMaximumIter', 'TDMaxSteps')
226 !%Variable TDMaxSteps
227 !%Type integer
228 !%Default 1500
229 !%Section Time-Dependent::Propagation
230 !%Description
231 !% Number of time-propagation steps that will be performed. You
232 !% cannot use this variable together with <tt>TDPropagationTime</tt>.
233 !%End
234 default = 1500
235 if (propagation_time > m_zero) default = nint(propagation_time/td%dt)
236 call parse_variable(namespace, 'TDMaxSteps', default, td%max_iter)
237
238 if (propagation_time <= m_zero) propagation_time = td%dt*td%max_iter
240 call messages_print_var_value('TDPropagationTime', propagation_time, unit = units_out%time, namespace=namespace)
241 call messages_print_var_value('TDMaxSteps', td%max_iter, namespace=namespace)
242
243 if (td%max_iter < 1) then
244 write(message(1), '(a,i6,a)') "Input: '", td%max_iter, "' is not a valid value for TDMaxSteps."
245 message(2) = '(TDMaxSteps <= 1)'
246 call messages_fatal(2, namespace=namespace)
247 end if
248
249 td%iter = 0
250
251 td%dt = td%dt
252
253 lasers => list_get_lasers(ext_partners)
254
255 ! now the photoelectron stuff
256 call pes_init(td%pesv, namespace, space, gr, gr%box, st, outp%restart_write_interval, hm%kpoints, &
257 hm%abs_boundaries, ext_partners, td%max_iter, td%dt)
258
259 !%Variable TDDynamics
260 !%Type integer
261 !%Default ehrenfest
262 !%Section Time-Dependent::Propagation
263 !%Description
264 !% Type of dynamics to follow during a time propagation.
265 !% For BO, you must set <tt>MoveIons = yes</tt>.
266 !%Option ehrenfest 1
267 !% Ehrenfest dynamics.
268 !%Option bo 2
269 !% Born-Oppenheimer (Experimental).
270 !%End
271
272 call parse_variable(namespace, 'TDDynamics', ehrenfest, td%dynamics)
273 if (.not. varinfo_valid_option('TDDynamics', td%dynamics)) call messages_input_error(namespace, 'TDDynamics')
274 call messages_print_var_option('TDDynamics', td%dynamics, namespace=namespace)
275 if (td%dynamics .ne. ehrenfest) then
276 if (.not. td%ions_dyn%is_active()) then
277 message(1) = "TDDynamics=bo can only be used if MoveIons=yes or CellDynamics=yes"
278 call messages_fatal(1, namespace=namespace)
279 end if
280 end if
281
282 !%Variable RecalculateGSDuringEvolution
283 !%Type logical
284 !%Default no
285 !%Section Time-Dependent::Propagation
286 !%Description
287 !% In order to calculate some information about the system along the
288 !% evolution (e.g. projection onto the ground-state KS determinant,
289 !% projection of the TDKS spin-orbitals onto the ground-state KS
290 !% spin-orbitals), the ground-state KS orbitals are needed. If the
291 !% ionic potential changes -- that is, the ions move -- one may want
292 !% to recalculate the ground state. You may do this by setting this
293 !% variable.
294 !%
295 !% The recalculation is not done every time step, but only every
296 !% <tt>RestartWriteInterval</tt> time steps.
297 !%End
298 call parse_variable(namespace, 'RecalculateGSDuringEvolution', .false., td%recalculate_gs)
299 if (hm%lda_u_level /= dft_u_none .and. td%recalculate_gs) then
300 call messages_not_implemented("DFT+U with RecalculateGSDuringEvolution=yes", namespace=namespace)
301 end if
302
303 !%Variable TDScissor
304 !%Type float
305 !%Default 0.0
306 !%Section Time-Dependent
307 !%Description
308 !% (experimental) If set, a scissor operator will be applied in the
309 !% Hamiltonian, shifting the excitation energies by the amount
310 !% specified. By default, it is not applied.
311 !%End
312 call parse_variable(namespace, 'TDScissor', m_zero, td%scissor)
313 call messages_print_var_value('TDScissor', td%scissor, namespace=namespace)
314
315 call propagator_elec_init(gr, namespace, st, td%tr, hm%ks_pot, td%ions_dyn%is_active() .and.&
316 list_has_gauge_field(ext_partners), family_is_mgga_with_exc(ks%xc), td%ions_dyn%cell_relax())
317
318 if (associated(lasers) .and. mpi_world%is_root()) then
319 call messages_print_with_emphasis(msg="Time-dependent external fields", namespace=namespace)
320 call laser_write_info(lasers%lasers, dt=td%dt, max_iter=td%max_iter, namespace=namespace)
321 call messages_print_with_emphasis(namespace=namespace)
322 end if
323
324 !%Variable TDEnergyUpdateIter
325 !%Type integer
326 !%Section Time-Dependent::Propagation
327 !%Description
328 !% This variable controls after how many iterations Octopus
329 !% updates the total energy during a time-propagation run. For
330 !% iterations where the energy is not updated, the last calculated
331 !% value is reported. If you set this variable to 1, the energy
332 !% will be calculated in each step.
333 !%End
334
335 default = 10
336 call parse_variable(namespace, 'TDEnergyUpdateIter', default, td%energy_update_iter)
337
338 if (gr%der%boundaries%spiralBC .and. hm%ep%reltype == spin_orbit) then
339 message(1) = "Generalized Bloch theorem cannot be used with spin-orbit coupling."
340 call messages_fatal(1, namespace=namespace)
341 end if
342
343 if (gr%der%boundaries%spiralBC) then
344 if (any(abs(hm%kick%easy_axis(1:2)) > m_epsilon)) then
345 message(1) = "Generalized Bloch theorem cannot be used for an easy axis not along the z direction."
346 call messages_fatal(1, namespace=namespace)
347 end if
348 end if
349
350 !%Variable TDFreezeOrbitals
351 !%Type integer
352 !%Default 0
353 !%Section Time-Dependent
354 !%Description
355 !% (Experimental) You have the possibility of "freezing" a number of orbitals during a time-propagation.
356 !% The Hartree and exchange-correlation potential due to these orbitals (which
357 !% will be the lowest-energy ones) will be added during the propagation, but the orbitals
358 !% will not be propagated.
359 !%Option sae -1
360 !% Single-active-electron approximation. This option is only valid for time-dependent
361 !% calculations (<tt>CalculationMode = td</tt>). Also, the nuclei should not move.
362 !% The idea is that all orbitals except the last one are frozen. The orbitals are to
363 !% be read from a previous ground-state calculation. The active orbital is then treated
364 !% as independent (whether it contains one electron or two) -- although it will
365 !% feel the Hartree and exchange-correlation potentials from the ground-state electronic
366 !% configuration.
367 !%
368 !% It is almost equivalent to setting <tt>TDFreezeOrbitals = N-1</tt>, where <tt>N</tt> is the number
369 !% of orbitals, but not completely.
370 !%End
371 call parse_variable(namespace, 'TDFreezeOrbitals', 0, td%freeze_orbitals)
372
373 if (td%freeze_orbitals /= 0) then
374 call messages_experimental('TDFreezeOrbitals', namespace=namespace)
375
376 if (hm%lda_u_level /= dft_u_none) then
377 call messages_not_implemented('TDFreezeOrbitals with DFT+U', namespace=namespace)
378 end if
379 end if
380
381 !%Variable TDDMPropagation
382 !%Type integer
383 !%Default no_propagation
384 !%Section Time-Dependent
385 !%Description
386 !% Decides whether to propagate the density matrix, in the spectral basis,
387 !% together with the KS orbitals.
388 !% The resulting time-depented orbilas are to be considered the natural orbitals
389 !% of the denstiy matrix and with the corresponding eigenvalues the occupation numbers.
390 !% In order to work properly this method requires an number of empty states
391 !% specified with the option "ExtraStates".
392 !% This is an experimental feature.
393 !%Option no_propagation 01
394 !% No density time propagation.
395 !%Option master_equation 02
396 !% Density time propagation with the Lindblad master equation.
397 !%End
398
399 call parse_variable(global_namespace, 'TDDMPropagation', option__tddmpropagation__no_propagation, td%dmp%calculation_mode)
400 if (.not. varinfo_valid_option('TDDMPropagation', td%dmp%calculation_mode)) then
401 call messages_input_error(global_namespace, 'TDDMPropagation')
402 endif
403 call messages_print_var_option('TDDMPropagation', td%dmp%calculation_mode, namespace=namespace)
404
405 if (td%dmp%calculation_mode /= option__tddmpropagation__no_propagation) then
406 call messages_experimental('TDDMPropagation', namespace=namespace)
407 call td%dmp%init(namespace, st)
408 end if
409
410 pop_sub(td_init)
411 nullify(lasers)
412
413 end subroutine td_init
414
415 ! ---------------------------------------------------------
416 subroutine td_init_run(td, namespace, mc, gr, ions, st, ks, hm, ext_partners, outp, space, from_scratch)
417 type(td_t), intent(inout) :: td
418 type(namespace_t), intent(in) :: namespace
419 type(multicomm_t), intent(inout) :: mc
420 type(grid_t), intent(inout) :: gr
421 type(ions_t), intent(inout) :: ions
422 type(states_elec_t), intent(inout) :: st
423 type(v_ks_t), intent(inout) :: ks
424 type(hamiltonian_elec_t), intent(inout) :: hm
425 type(partner_list_t), intent(in) :: ext_partners
426 type(output_t), intent(inout) :: outp
427 type(electron_space_t), intent(in) :: space
428 logical, intent(inout) :: from_scratch
429 push_sub(td_init_run)
430
431 ! NOTE: please do not change code in this function, but only in functions
432 ! called from here because the logic of this function is replicated in the
433 ! multisystem framework in different places
434
435 call td_allocate_wavefunctions(td, namespace, mc, gr, ions, st, hm, space)
436 call td_init_gaugefield(td, namespace, gr, st, ks, hm, ext_partners, space)
437
438 td%from_scratch = from_scratch
439
440 if (.not. td%from_scratch) then
441 call td_load_restart_from_td(td, namespace, space, mc, gr, ext_partners, st, ks, hm, td%from_scratch)
442 if (td%from_scratch) then
443 message(1) = "Unable to read time-dependent restart information: Starting from scratch"
444 call messages_warning(1, namespace=namespace)
445 end if
446 end if
447
448 if (td%iter >= td%max_iter) then
449 message(1) = "All requested iterations have already been done. Use FromScratch = yes if you want to redo them."
450 call messages_info(1, namespace=namespace)
452 td%iter = td%iter + 1
453 if ((td%ions_dyn%is_active()) .and. td%recalculate_gs) call td%restart_load%end()
454 pop_sub(td_init_run)
455 return
456 end if
457
458 if (td%from_scratch) then
459 call td_load_restart_from_gs(td, namespace, space, mc, gr, ext_partners, st, ks, hm)
460 end if
461
462 call td_init_with_wavefunctions(td, namespace, space, mc, gr, ions, ext_partners, st, ks, hm, outp, td%from_scratch)
463
464 if (td%dmp%calculation_mode /= option__tddmpropagation__no_propagation) then
465 call dm_propagation_init_run(td%dmp%adiabatic_st, namespace, space, gr, st, hm, mc)
466 end if
467
468 pop_sub(td_init_run)
469 end subroutine td_init_run
470
471 ! ---------------------------------------------------------
472 subroutine td_allocate_wavefunctions(td, namespace, mc, gr, ions, st, hm, space)
473 type(td_t), intent(inout) :: td
474 type(namespace_t), intent(in) :: namespace
475 type(multicomm_t), intent(inout) :: mc
476 type(grid_t), intent(inout) :: gr
477 type(ions_t), intent(inout) :: ions
478 type(states_elec_t), intent(inout) :: st
479 type(hamiltonian_elec_t), intent(inout) :: hm
480 class(space_t), intent(in) :: space
481
483
484 ! Allocate wavefunctions during time-propagation
485 if (td%dynamics == ehrenfest) then
486 !Note: this is not really clean to do this
487 if (hm%lda_u_level /= dft_u_none .and. states_are_real(st)) then
488 call lda_u_end(hm%lda_u)
489 !complex wfs are required for Ehrenfest
490 call states_elec_allocate_wfns(st, gr, type_cmplx, packed=.true.)
491 call lda_u_init(hm%lda_u, namespace, space, hm%lda_u_level, gr, ions, st, mc, hm%kpoints)
492 else
493 !complex wfs are required for Ehrenfest
494 call states_elec_allocate_wfns(st, gr, type_cmplx, packed=.true.)
495 end if
496 else
497 call states_elec_allocate_wfns(st, gr, packed=.true.)
498 call scf_init(td%scf, namespace, gr, ions, st, mc, hm, space)
499 end if
500
502 end subroutine td_allocate_wavefunctions
503
504 ! ---------------------------------------------------------
505 subroutine td_init_gaugefield(td, namespace, gr, st, ks, hm, ext_partners, space)
506 type(td_t), intent(inout) :: td
507 type(namespace_t), intent(in) :: namespace
508 type(grid_t), intent(inout) :: gr
509 type(states_elec_t), intent(inout) :: st
510 type(v_ks_t), intent(inout) :: ks
511 type(hamiltonian_elec_t), intent(inout) :: hm
512 type(partner_list_t), intent(in) :: ext_partners
513 class(space_t), intent(in) :: space
514
515 type(gauge_field_t), pointer :: gfield
516
517 push_sub(td_init_gaugefield)
518
519 gfield => list_get_gauge_field(ext_partners)
520 if(associated(gfield)) then
521 if (gauge_field_is_used(gfield)) then
522 !if the gauge field is applied, we need to tell v_ks to calculate the current
523 call v_ks_calculate_current(ks, .true.)
524
525 ! initialize the vector field and update the hamiltonian
526 call gauge_field_init_vec_pot(gfield, st%qtot)
527 call hm%update(gr, namespace, space, ext_partners, time = td%dt*td%iter)
528
529 end if
530 end if
531
532 pop_sub(td_init_gaugefield)
533 end subroutine td_init_gaugefield
534
535 ! ---------------------------------------------------------
536 subroutine td_end(td)
537 type(td_t), intent(inout) :: td
538
539 push_sub(td_end)
540
541 call pes_end(td%pesv)
542 call propagator_elec_end(td%tr) ! clean the evolution method
543 call ion_dynamics_end(td%ions_dyn)
544
545 if (td%dmp%calculation_mode /= option__tddmpropagation__no_propagation) then
546 call states_elec_end(td%dmp%adiabatic_st)
547 endif
548
549 if (td%dynamics == bo) call scf_end(td%scf)
550 pop_sub(td_end)
551
552 end subroutine td_end
553
554 ! ---------------------------------------------------------
555 subroutine td_end_run(td, st, hm)
556 type(td_t), intent(inout) :: td
557 type(states_elec_t), intent(inout) :: st
558 type(hamiltonian_elec_t), intent(inout) :: hm
559
560 push_sub(td_end_run)
561
562 if (st%pack_states .and. hm%apply_packed()) call st%unpack()
563
564 call td%restart_dump%end()
565 call td_write_end(td%write_handler)
566
567 ! free memory
569 if ((td%ions_dyn%is_active()).and. td%recalculate_gs) call td%restart_load%end()
570
571 pop_sub(td_end_run)
572 end subroutine td_end_run
573
574 ! ---------------------------------------------------------
575 subroutine td_run(td, namespace, mc, gr, ions, st, ks, hm, ext_partners, outp, space, from_scratch)
576 type(td_t), intent(inout) :: td
577 type(namespace_t), intent(in) :: namespace
578 type(multicomm_t), intent(inout) :: mc
579 type(grid_t), intent(inout) :: gr
580 type(ions_t), intent(inout) :: ions
581 type(states_elec_t), intent(inout) :: st
582 type(v_ks_t), intent(inout) :: ks
583 type(hamiltonian_elec_t), intent(inout) :: hm
584 type(partner_list_t), intent(in) :: ext_partners
585 type(output_t), intent(inout) :: outp
586 type(electron_space_t), intent(in) :: space
587 logical, intent(inout) :: from_scratch
588
589 logical :: stopping
590 integer :: iter, scsteps
591 real(real64) :: etime
592 real(real64) :: wall_time, simulation_time, speed_fs_per_day
593 character(len=20) :: fmt
594
595 push_sub(td_run)
596
597 etime = loct_clock()
598 ! This is the time-propagation loop. It starts at t=0 and finishes
599 ! at td%max_iter*dt. The index i runs from 1 to td%max_iter, and
600 ! step "iter" means propagation from (iter-1)*dt to iter*dt.
601 propagation: do iter = td%iter, td%max_iter
602
603 stopping = clean_stop(mc%master_comm) .or. walltimer_alarm(mc%master_comm)
604
605 call profiling_in("TIME_STEP")
606
607 if (iter > 1) then
608 if (((iter-1)*td%dt <= hm%kick%time) .and. (iter*td%dt > hm%kick%time)) then
609 if (.not. hm%pcm%localf) then
610 call kick_apply(space, gr, st, td%ions_dyn, ions, hm%kick, hm%psolver, hm%kpoints)
611 else
612 call kick_apply(space, gr, st, td%ions_dyn, ions, hm%kick, hm%psolver, hm%kpoints, pcm = hm%pcm)
613 end if
614 call td_write_kick(outp, namespace, space, gr, hm%kick, ions, iter)
615 !We activate the sprial BC only after the kick,
616 !to be sure that the first iteration corresponds to the ground state
617 if (gr%der%boundaries%spiralBC) gr%der%boundaries%spiral = .true.
618 end if
619 end if
620
621 ! time iterate the system, one time step.
622 select case (td%dynamics)
623 case (ehrenfest)
624 call propagator_elec_dt(ks, namespace, space, hm, gr, st, td%tr, iter*td%dt, td%dt, iter, td%ions_dyn, &
625 ions, ext_partners, mc, outp, td%write_handler, scsteps = scsteps, &
626 update_energy = (mod(iter, td%energy_update_iter) == 0) .or. (iter == td%max_iter))
627 case (bo)
628 call propagator_elec_dt_bo(td%scf, namespace, space, gr, ks, st, hm, ions, ext_partners, mc, iter, &
629 td%dt, td%ions_dyn, scsteps)
630 end select
632 if (td%dmp%calculation_mode /= option__tddmpropagation__no_propagation) then
633 call dm_propagation_run(td%dmp, namespace, space, gr, ions, st, mc, hm, ks, iter, td%dt, ext_partners, &
634 update_energy = (mod(iter, td%energy_update_iter) == 0) .or. (iter == td%max_iter))
635 end if
636
637 !Apply mask absorbing boundaries
638 if (hm%abs_boundaries%abtype == mask_absorbing) then
639 if (states_are_real(st)) then
640 call dvmask(gr, hm, st)
641 else
642 call zvmask(gr, hm, st)
643 end if
644 end if
645
646 !Photoelectron stuff
647 if (td%pesv%calc_spm .or. td%pesv%calc_mask .or. td%pesv%calc_flux) then
648 call pes_calc(td%pesv, namespace, space, gr, st, td%dt, iter, gr%der, hm%kpoints, ext_partners, stopping)
649 end if
651 call td_write_iter(td%write_handler, namespace, space, outp, gr, st, hm, ions, ext_partners, &
652 hm%kick, ks, td%dt, iter, mc, td%recalculate_gs)
653
654 ! write down data
655 call td_check_point(td, namespace, mc, gr, ions, st, ks, hm, ext_partners, outp, space, &
656 iter, scsteps, etime, stopping, from_scratch)
657
658 ! check if debug mode should be enabled or disabled on the fly
659 call io_debug_on_the_fly(namespace)
660
661 call profiling_out("TIME_STEP")
662 if (stopping) exit
663
664 end do propagation
665
666 ! Output propagation speed in fs/day
667 wall_time = loct_clock() - walltimer_get_start_time()
668 simulation_time = td%dt * (iter - td%iter + 1)
669 speed_fs_per_day = units_from_atomic(unit_femtosecond, simulation_time) / (wall_time / 86400.0_real64)
670 if (speed_fs_per_day > 1e4 .or. speed_fs_per_day < 1e-3) then
671 fmt = "ES10.3"
672 else
673 fmt = "F10.3"
674 end if
675 write(message(1), '(a,'//trim(fmt)//',a)') 'Propagation speed: ', speed_fs_per_day, ' fs/day'
676 write(message(2), '(a)') ''
677 call messages_info(2, namespace=namespace)
678
679 pop_sub(td_run)
680 end subroutine td_run
681
682 subroutine td_print_header(namespace)
683 type(namespace_t), intent(in) :: namespace
684
685 push_sub(td_print_header)
686
687 write(message(1), '(a7,1x,a14,a14,a10,a17)') 'Iter ', 'Time ', 'Energy ', 'SC Steps', 'Elapsed Time '
688
689 call messages_info(1, namespace=namespace)
690 call messages_print_with_emphasis(namespace=namespace)
691
692 pop_sub(td_print_header)
693 end subroutine td_print_header
694
695 ! ---------------------------------------------------------
696 subroutine td_check_point(td, namespace, mc, gr, ions, st, ks, hm, ext_partners, outp, space, &
697 iter, scsteps, etime, stopping, from_scratch)
698 type(td_t), intent(inout) :: td
699 type(namespace_t), intent(in) :: namespace
700 type(multicomm_t), intent(in) :: mc
701 type(grid_t), intent(inout) :: gr
702 type(ions_t), intent(inout) :: ions
703 type(states_elec_t), intent(inout) :: st
704 type(v_ks_t), intent(inout) :: ks
705 type(hamiltonian_elec_t), intent(inout) :: hm
706 type(partner_list_t), intent(in) :: ext_partners
707 type(output_t), intent(in) :: outp
708 type(electron_space_t), intent(in) :: space
709 integer, intent(in) :: iter
710 integer, intent(in) :: scsteps
711 real(real64), intent(inout) :: etime
712 logical, intent(in) :: stopping
713 logical, intent(inout) :: from_scratch
714
715 integer :: ierr
716
717 push_sub(td_check_point)
718
719 call td_print_message(td, namespace, ions, hm, iter, scsteps, etime)
720
721 if (outp%anything_now(iter)) then ! output
722 call td_write_output(namespace, space, gr, st, hm, ks, outp, ions, ext_partners, iter, td%dt)
723 end if
724
725 if (mod(iter, outp%restart_write_interval) == 0 .or. iter == td%max_iter .or. stopping) then ! restart
726 !if (iter == td%max_iter) outp%iter = ii - 1
727 call td_write_data(td%write_handler)
728 call td_dump(td, namespace, space, gr, st, hm, ks, ext_partners, iter, ierr)
729 if (ierr /= 0) then
730 message(1) = "Unable to write time-dependent restart information."
731 call messages_warning(1, namespace=namespace)
732 end if
733
734 call pes_output(td%pesv, namespace, space, gr, st, iter, outp, td%dt, ions)
735
736 if ((td%ions_dyn%is_active()) .and. td%recalculate_gs) then
737 call messages_print_with_emphasis(msg='Recalculating the ground state.', namespace=namespace)
738 from_scratch = .false.
740 call electrons_ground_state_run(namespace, mc, gr, ions, ext_partners, st, ks, hm, outp, space, from_scratch)
741 call states_elec_allocate_wfns(st, gr, packed=.true.)
742 call td_load(td%restart_load, namespace, space, gr, st, hm, ext_partners, td, ks, ierr)
743 if (ierr /= 0) then
744 message(1) = "Unable to load TD states."
745 call messages_fatal(1, namespace=namespace)
746 end if
747 call density_calc(st, gr, st%rho)
748 call v_ks_calc(ks, namespace, space, hm, st, ions, ext_partners, &
749 calc_eigenval=.true., time = iter*td%dt, calc_energy=.true.)
750 call forces_calculate(gr, namespace, ions, hm, ext_partners, st, ks, t = iter*td%dt, dt = td%dt)
751 assert(.not. td%ions_dyn%cell_relax())
752 call messages_print_with_emphasis(msg="Time-dependent simulation proceeds", namespace=namespace)
753 call td_print_header(namespace)
754 end if
755 end if
756
757 pop_sub(td_check_point)
758 end subroutine td_check_point
759
760 ! ---------------------------------------------------------
761 subroutine td_print_message(td, namespace, ions, hm, iter, scsteps, etime)
762 type(td_t), intent(inout) :: td
763 type(namespace_t), intent(in) :: namespace
764 type(ions_t), intent(inout) :: ions
765 type(hamiltonian_elec_t), intent(inout) :: hm
766 integer, intent(in) :: iter
767 integer, intent(in) :: scsteps
768 real(real64), intent(inout) :: etime
769
770 push_sub(td_print_message)
771
772 write(message(1), '(i7,1x,2f14.6,i10,f14.3)') iter, units_from_atomic(units_out%time, iter*td%dt), &
773 units_from_atomic(units_out%energy, hm%energy%total + ions%kinetic_energy), scsteps, &
774 loct_clock() - etime
775 call messages_info(1, namespace=namespace)
776 call td_update_elapsed_time(etime)
778 pop_sub(td_print_message)
779 end subroutine td_print_message
780
781 ! ---------------------------------------------------------
782 subroutine td_update_elapsed_time(etime)
783 real(real64), intent(inout) :: etime
784
785 push_sub(td_update_elapsed_time)
786
787 etime = loct_clock()
788
790 end subroutine td_update_elapsed_time
792 ! ---------------------------------------------------------
793 subroutine td_init_with_wavefunctions(td, namespace, space, mc, gr, ions, ext_partners, st, ks, hm, outp, from_scratch)
794 type(td_t), intent(inout) :: td
795 type(namespace_t), intent(in) :: namespace
796 type(electron_space_t), intent(in) :: space
797 type(multicomm_t), intent(in) :: mc
798 type(grid_t), intent(inout) :: gr
799 type(ions_t), intent(inout) :: ions
800 type(partner_list_t), intent(in) :: ext_partners
801 type(states_elec_t), target, intent(inout) :: st
802 type(v_ks_t), intent(inout) :: ks
803 type(hamiltonian_elec_t), intent(inout) :: hm
804 type(output_t), intent(inout) :: outp
805 logical, intent(in) :: from_scratch
806
807 integer :: ierr
808 real(real64) :: x
809 real(real64) :: ndinitial(space%dim)
810 logical :: freeze_hxc, freeze_occ, freeze_u
811 type(restart_t) :: restart, restart_frozen
812 type(gauge_field_t), pointer :: gfield
813 type(lasers_t), pointer :: lasers
815
816 !We activate the sprial BC only after the kick,
817 !to be sure that the first iteration corresponds to the ground state
818 if (gr%der%boundaries%spiralBC) then
819 if ((td%iter-1)*td%dt > hm%kick%time) then
820 gr%der%boundaries%spiral = .true.
821 end if
822 hm%vnl%spin => st%spin
823 hm%phase%spin => st%spin
824 !We fill st%spin. In case of restart, we read it in td_load
825 if (from_scratch) call states_elec_fermi(st, namespace, gr)
826 end if
827
828 if (from_scratch) then
829 ! Initialize the occupation matrices and U for DFT+U
830 ! This must be called before parsing TDFreezeOccupations and TDFreezeU
831 ! in order that the code does properly the initialization.
832 call lda_u_update_occ_matrices(hm%lda_u, namespace, gr, st, hm%phase, hm%energy)
833 end if
834
835 if (td%freeze_orbitals > 0) then
836 if (from_scratch) then
837 ! In this case, we first freeze the orbitals, then calculate the Hxc potential.
838 call states_elec_freeze_orbitals(st, namespace, space, gr, mc, hm%kpoints, &
839 td%freeze_orbitals, family_is_mgga(ks%xc_family))
840 else
841 call restart%init(namespace, restart_td, restart_type_load, mc, ierr, mesh=gr)
842 if (ierr == 0) then
843 call td_load_frozen(namespace, restart, space, gr, st, hm, ierr)
844 end if
845 if (ierr /= 0) then
846 td%iter = 0
847 message(1) = "Unable to read frozen restart information."
848 call messages_fatal(1, namespace=namespace)
849 end if
850 call restart%end()
851 end if
852 write(message(1),'(a,i4,a,i4,a)') 'Info: The lowest', td%freeze_orbitals, &
853 ' orbitals have been frozen.', st%nst, ' will be propagated.'
854 call messages_info(1, namespace=namespace)
856 call density_calc(st, gr, st%rho)
857 call v_ks_calc(ks, namespace, space, hm, st, ions, ext_partners, calc_eigenval=.true., time = td%iter*td%dt)
858 else if (td%freeze_orbitals < 0) then
859 ! This means SAE approximation. We calculate the Hxc first, then freeze all
860 ! orbitals minus one.
861 write(message(1),'(a)') 'Info: The single-active-electron approximation will be used.'
862 call messages_info(1, namespace=namespace)
863 call v_ks_calc(ks, namespace, space, hm, st, ions, ext_partners, calc_eigenval=.true., time = td%iter*td%dt)
864 if (from_scratch) then
865 call states_elec_freeze_orbitals(st, namespace, space, gr, mc, hm%kpoints, st%nst-1, family_is_mgga(ks%xc_family))
866 else
867 call messages_not_implemented("TDFreezeOrbials < 0 with FromScratch=no", namespace=namespace)
868 end if
870 call v_ks_freeze_hxc(ks)
871 call density_calc(st, gr, st%rho)
872 else
873 ! Normal run.
874 call density_calc(st, gr, st%rho)
875 call v_ks_calc(ks, namespace, space, hm, st, ions, ext_partners, calc_eigenval=.true., time = td%iter*td%dt)
876 end if
878 !%Variable TDFreezeHXC
879 !%Type logical
880 !%Default no
881 !%Section Time-Dependent
882 !%Description
883 !% The electrons are evolved as independent particles feeling the Hartree and
884 !% exchange-correlation potentials from the ground-state electronic configuration.
885 !%End
886 call parse_variable(namespace, 'TDFreezeHXC', .false., freeze_hxc)
887 if (freeze_hxc) then
888 write(message(1),'(a)') 'Info: Freezing Hartree and exchange-correlation potentials.'
889 call messages_info(1, namespace=namespace)
890
891 if (.not. from_scratch) then
892
893 call restart_frozen%init(namespace, restart_gs, restart_type_load, mc, ierr, mesh=gr, exact=.true.)
894 call states_elec_load(restart_frozen, namespace, space, st, gr, hm%kpoints, ierr, label = ": gs")
895 call states_elec_transform(st, namespace, space, restart_frozen, gr, hm%kpoints)
896 call restart_frozen%end()
897
898 call density_calc(st, gr, st%rho)
899 call v_ks_calc(ks, namespace, space, hm, st, ions, ext_partners, calc_eigenval=.true., time = td%iter*td%dt)
900
901 call restart_frozen%init(namespace, restart_td, restart_type_load, mc, ierr, mesh=gr)
902 call states_elec_load(restart_frozen, namespace, space, st, gr, hm%kpoints, ierr, iter=td%iter, label = ": td")
903 call restart_frozen%end()
904 call hm%ks_pot%run_zero_iter(td%tr%vks_old)
905
906 end if
907
908 call v_ks_freeze_hxc(ks)
909
910 end if
911
912 x = minval(st%eigenval(st%st_start, :))
913 if (st%parallel_in_states) then
914 call st%mpi_grp%bcast(x, 1, mpi_double_precision, 0)
915 end if
916 call hm%update_span(gr%spacing(1:space%dim), x, namespace)
917 ! initialize Fermi energy
918 call states_elec_fermi(st, namespace, gr, compute_spin = .not. gr%der%boundaries%spiralBC)
919 call energy_calc_total(namespace, space, hm, gr, st, ext_partners)
920
921 !%Variable TDFreezeDFTUOccupations
922 !%Type logical
923 !%Default no
924 !%Section Time-Dependent
925 !%Description
926 !% The occupation matrices than enters in the DFT+U potential
927 !% are not evolved during the time evolution.
928 !%End
929 call parse_variable(namespace, 'TDFreezeDFTUOccupations', .false., freeze_occ)
930 if (freeze_occ) then
931 write(message(1),'(a)') 'Info: Freezing DFT+U occupation matrices that enters in the DFT+U potential.'
932 call messages_info(1, namespace=namespace)
933 call lda_u_freeze_occ(hm%lda_u)
934
935 !In this case we should reload GS wavefunctions
936 if (hm%lda_u_level /= dft_u_none .and. .not. from_scratch) then
937 call restart_frozen%init(namespace, restart_gs, restart_type_load, mc, ierr, mesh=gr)
938 call lda_u_load(restart_frozen, hm%lda_u, st, hm%energy%dft_u, ierr, occ_only = .true.)
939 call restart_frozen%end()
940 end if
941 end if
942
943 !%Variable TDFreezeU
944 !%Type logical
945 !%Default no
946 !%Section Time-Dependent
947 !%Description
948 !% The effective U of DFT+U is not evolved during the time evolution.
949 !%End
950 call parse_variable(namespace, 'TDFreezeU', .false., freeze_u)
951 if (freeze_u) then
952 write(message(1),'(a)') 'Info: Freezing the effective U of DFT+U.'
953 call messages_info(1, namespace=namespace)
954 call lda_u_freeze_u(hm%lda_u)
955
956 !In this case we should reload GS wavefunctions
957 if (hm%lda_u_level == dft_u_acbn0 .and. .not. from_scratch) then
958 call restart_frozen%init(namespace, restart_gs, restart_type_load, mc, ierr, mesh=gr)
959 call lda_u_load(restart_frozen, hm%lda_u, st, hm%energy%dft_u, ierr, u_only = .true.)
960 call restart_frozen%end()
961 write(message(1),'(a)') 'Loaded GS effective U of DFT+U'
962 call messages_info(1, namespace=namespace)
963 call lda_u_write_u(hm%lda_u, namespace=namespace)
964 call lda_u_write_v(hm%lda_u, namespace=namespace)
965 end if
966 end if
967
968 ! This needs to be called before the calculation of the forces,
969 ! as we need to test of we output the forces or not
970 call td_write_init(td%write_handler, namespace, space, outp, gr, st, hm, ions, ext_partners, &
971 ks, td%ions_dyn%is_active(), &
972 list_has_gauge_field(ext_partners), hm%kick, td%iter, td%max_iter, td%dt, mc)
973
974 ! Resets the nondipole integration after laser-file has been written.
975 lasers => list_get_lasers(ext_partners)
976 if(associated(lasers)) then
977 if (lasers_with_nondipole_field(lasers)) then
978 ndinitial(1:space%dim)=m_zero
979 call lasers_set_nondipole_parameters(lasers,ndinitial,m_zero)
980 end if
981 end if
982 nullify(lasers)
983
984 call td_init_ions_and_forces(td, namespace, space, gr, ions, ext_partners, st, ks, hm, outp)
985
986 if (td%scissor > m_epsilon) then
987 call scissor_init(hm%scissor, namespace, space, st, gr, hm%kpoints, hm%phase, td%scissor, mc)
988 end if
989
990 if (td%iter == 0) call td_run_zero_iter(td, namespace, space, gr, ions, st, ks, hm, ext_partners, outp, mc)
991
992 gfield => list_get_gauge_field(ext_partners)
993 if(associated(gfield)) then
994 if (gauge_field_is_propagated(gfield)) then
995 call gauge_field_get_force(gfield, gr, st%d%spin_channels, st%current, ks%xc%lrc)
996 if (abs(ks%xc%lrc%alpha) > m_epsilon) then
997 call messages_experimental('TD-LRC kernel')
998 end if
999 endif
1000 end if
1001
1002 !call td_check_trotter(td, sys, h)
1003 td%iter = td%iter + 1
1004
1005 call td%restart_dump%init(namespace, restart_td, restart_type_dump, mc, ierr, mesh=gr)
1006 if (td%ions_dyn%is_active() .and. td%recalculate_gs) then
1007 ! We will also use the TD restart directory as temporary storage during the time propagation
1008 call td%restart_load%init(namespace, restart_td, restart_type_load, mc, ierr, mesh=gr)
1009 end if
1010
1011 call messages_print_with_emphasis(msg="Time-Dependent Simulation", namespace=namespace)
1012 call td_print_header(namespace)
1013
1014 if (td%pesv%calc_spm .or. td%pesv%calc_mask .and. from_scratch) then
1015 call pes_init_write(td%pesv,gr,st, namespace)
1016 end if
1017
1018 if (st%pack_states .and. hm%apply_packed()) call st%pack()
1019
1021 end subroutine td_init_with_wavefunctions
1022
1023 ! ---------------------------------------------------------
1024 subroutine td_init_ions_and_forces(td, namespace, space, gr, ions, ext_partners, st, ks, hm, outp)
1025 type(td_t), intent(inout) :: td
1026 type(namespace_t), intent(in) :: namespace
1027 type(electron_space_t), intent(in) :: space
1028 type(grid_t), intent(inout) :: gr
1029 type(ions_t), intent(inout) :: ions
1030 type(partner_list_t), intent(in) :: ext_partners
1031 type(states_elec_t), target, intent(inout) :: st
1032 type(v_ks_t), intent(inout) :: ks
1033 type(hamiltonian_elec_t), intent(inout) :: hm
1034 type(output_t), intent(inout) :: outp
1035
1036 push_sub(td_init_ions_and_forces)
1037
1038 ! Calculate initial forces and kinetic energy
1039 if (td%ions_dyn%ions_move()) then
1040 if (td%iter > 0) then
1041 call td_read_coordinates(td, namespace, ions)
1042 if (ion_dynamics_drive_ions(td%ions_dyn)) then
1043 call ion_dynamics_propagate(td%ions_dyn, ions, td%iter*td%dt, td%dt, namespace)
1044 end if
1045 call hamiltonian_elec_epot_generate(hm, namespace, space, gr, ions, ext_partners, st, time = td%iter*td%dt)
1046 ! recompute potential because the ions have moved
1047 call v_ks_calc(ks, namespace, space, hm, st, ions, ext_partners, calc_eigenval=.true., time = td%iter*td%dt)
1048 end if
1049
1050 call forces_calculate(gr, namespace, ions, hm, ext_partners, st, ks, t = td%iter*td%dt, dt = td%dt)
1051
1052 call ions%update_kinetic_energy()
1053 else
1054 if (outp%what(option__output__forces) .or. td%write_handler%out(out_separate_forces)%write) then
1055 call forces_calculate(gr, namespace, ions, hm, ext_partners, st, ks, t = td%iter*td%dt, dt = td%dt)
1056 end if
1057 end if
1058
1059 if (outp%what(option__output__stress) .or. td%ions_dyn%cell_relax()) then
1060 call stress_calculate(namespace, gr, hm, st, ions, ks, ext_partners)
1061 if (td%ions_dyn%cell_relax()) then
1062 call td%ions_dyn%update_stress(ions%space, st%stress_tensors%total, ions%latt%rlattice, ions%latt%rcell_volume)
1063 end if
1064 end if
1065
1067 end subroutine td_init_ions_and_forces
1068
1069 ! ---------------------------------------------------------
1070 subroutine td_load_restart_from_td(td, namespace, space, mc, gr, ext_partners, st, ks, hm, from_scratch)
1071 type(td_t), intent(inout) :: td
1072 type(namespace_t), intent(in) :: namespace
1073 class(space_t), intent(in) :: space
1074 type(multicomm_t), intent(in) :: mc
1075 type(grid_t), intent(inout) :: gr
1076 type(partner_list_t), intent(in) :: ext_partners
1077 type(states_elec_t), target, intent(inout) :: st
1078 type(v_ks_t), intent(inout) :: ks
1079 type(hamiltonian_elec_t), intent(inout) :: hm
1080 logical, intent(inout) :: from_scratch
1081
1082 integer :: ierr
1083 type(restart_t) :: restart
1084
1085 push_sub(td_load_restart_from_td)
1086
1087 !We redistribute the states before the restarting
1088 if (td%freeze_orbitals > 0) then
1089 call states_elec_freeze_redistribute_states(st, namespace, gr, mc, td%freeze_orbitals)
1090 end if
1091
1092 call restart%init(namespace, restart_td, restart_type_load, mc, ierr, mesh=gr)
1093 if (ierr == 0) then
1094 call td_load(restart, namespace, space, gr, st, hm, ext_partners, td, ks, ierr)
1095 end if
1096 call restart%end()
1097 if (ierr /= 0) then
1098 from_scratch = .true.
1099 td%iter = 0
1100 end if
1101
1103 end subroutine td_load_restart_from_td
1104
1105 ! ---------------------------------------------------------
1106 subroutine td_load_restart_from_gs(td, namespace, space, mc, gr, ext_partners, st, ks, hm)
1107 type(td_t), intent(inout) :: td
1108 type(namespace_t), intent(in) :: namespace
1109 class(space_t), intent(in) :: space
1110 type(multicomm_t), intent(in) :: mc
1111 type(grid_t), intent(inout) :: gr
1112 type(partner_list_t), intent(in) :: ext_partners
1113 type(states_elec_t), target, intent(inout) :: st
1114 type(v_ks_t), intent(inout) :: ks
1115 type(hamiltonian_elec_t), intent(inout) :: hm
1116
1117 integer :: ierr
1118 type(restart_t) :: restart
1120 push_sub(td_load_restart_from_gs)
1121
1122 call restart%init(namespace, restart_gs, restart_type_load, mc, ierr, mesh=gr, exact=.true.)
1123
1124 if (.not. st%only_userdef_istates) then
1125 if (ierr == 0) then
1126 call states_elec_load(restart, namespace, space, st, gr, hm%kpoints, ierr, label = ": gs")
1127 end if
1128 if (ierr /= 0) then
1129 message(1) = 'Unable to read ground-state wavefunctions.'
1130 call messages_fatal(1, namespace=namespace)
1131 end if
1132 end if
1133
1134 ! check if we should deploy user-defined wavefunctions.
1135 ! according to the settings in the input file the routine
1136 ! overwrites orbitals that were read from restart/gs
1137 if (parse_is_defined(namespace, 'UserDefinedStates')) then
1138 call states_elec_read_user_def_orbitals(gr, namespace, space, st)
1139 end if
1140
1141 call states_elec_transform(st, namespace, space, restart, gr, hm%kpoints)
1142 call restart%end()
1143
1145 end subroutine td_load_restart_from_gs
1146
1147 ! ---------------------------------------------------------
1148 subroutine td_run_zero_iter(td, namespace, space, gr, ions, st, ks, hm, ext_partners, outp, mc)
1149 type(td_t), intent(inout) :: td
1150 type(namespace_t), intent(in) :: namespace
1151 type(electron_space_t), intent(in) :: space
1152 type(grid_t), intent(inout) :: gr
1153 type(ions_t), intent(inout) :: ions
1154 type(states_elec_t), intent(inout) :: st
1155 type(v_ks_t), intent(inout) :: ks
1156 type(hamiltonian_elec_t), intent(inout) :: hm
1157 type(partner_list_t), intent(in) :: ext_partners
1158 type(output_t), intent(in) :: outp
1159 type(multicomm_t), intent(in) :: mc
1160
1161 push_sub(td_run_zero_iter)
1162
1163 call td_write_iter(td%write_handler, namespace, space, outp, gr, st, hm, ions, ext_partners, &
1164 hm%kick, ks, td%dt, 0, mc, td%recalculate_gs)
1166 ! I apply the delta electric field *after* td_write_iter, otherwise the
1167 ! dipole matrix elements in write_proj are wrong
1168 if (abs(hm%kick%time) <= m_epsilon) then
1169 if (.not. hm%pcm%localf) then
1170 call kick_apply(space, gr, st, td%ions_dyn, ions, hm%kick, hm%psolver, hm%kpoints)
1171 else
1172 call kick_apply(space, gr, st, td%ions_dyn, ions, hm%kick, hm%psolver, hm%kpoints, pcm = hm%pcm)
1173 end if
1174 call td_write_kick(outp, namespace, space, gr, hm%kick, ions, 0)
1175
1176 !We activate the sprial BC only after the kick
1177 if (gr%der%boundaries%spiralBC) then
1178 gr%der%boundaries%spiral = .true.
1179 end if
1180 end if
1181 call hm%ks_pot%run_zero_iter(td%tr%vks_old)
1182
1183 if (any(outp%output_interval > 0)) then
1184 call td_write_data(td%write_handler)
1185 call td_write_output(namespace, space, gr, st, hm, ks, outp, ions, ext_partners, 0)
1186 end if
1187
1188 pop_sub(td_run_zero_iter)
1189 end subroutine td_run_zero_iter
1190
1191
1192 ! ---------------------------------------------------------
1194 subroutine td_read_coordinates(td, namespace, ions)
1195 type(td_t), intent(in) :: td
1196 type(namespace_t), intent(in) :: namespace
1197 type(ions_t), intent(inout) :: ions
1198
1199 integer :: iatom, iter, iunit
1200
1202
1203 iunit = io_open('td.general/coordinates', namespace, action='read', status='old', die=.false.)
1204 if (iunit == -1) then
1205 message(1) = "Could not open file '"//trim(io_workpath('td.general/coordinates', namespace))//"'."
1206 message(2) = "Starting simulation from initial geometry."
1207 call messages_warning(2, namespace=namespace)
1208 pop_sub(td_read_coordinates)
1209 return
1210 end if
1211
1212 call io_skip_header(iunit)
1213 do iter = 0, td%iter - 1
1214 read(iunit, *) ! skip previous iterations... sorry, but no portable seek in Fortran
1215 end do
1216 read(iunit, '(32x)', advance='no') ! skip the time index.
1217
1218 do iatom = 1, ions%natoms
1219 read(iunit, '(3es24.16)', advance='no') ions%pos(:, iatom)
1220 ions%pos(:, iatom) = units_to_atomic(units_out%length, ions%pos(:, iatom))
1221 end do
1222 do iatom = 1, ions%natoms
1223 read(iunit, '(3es24.16)', advance='no') ions%vel(:, iatom)
1224 ions%vel(:, iatom) = units_to_atomic(units_out%velocity, ions%vel(:, iatom))
1225 end do
1226 do iatom = 1, ions%natoms
1227 read(iunit, '(3es24.16)', advance='no') ions%tot_force(:, iatom)
1228 ions%tot_force(:, iatom) = units_to_atomic(units_out%force, ions%tot_force(:, iatom))
1229 end do
1230
1231 call io_close(iunit)
1232
1233 pop_sub(td_read_coordinates)
1234 end subroutine td_read_coordinates
1235
1236 ! ---------------------------------------------------------
1237 subroutine td_dump(td, namespace, space, gr, st, hm, ks, ext_partners, iter, ierr)
1238 type(td_t), intent(in) :: td
1239 type(namespace_t), intent(in) :: namespace
1240 class(space_t), intent(in) :: space
1241 type(grid_t), intent(in) :: gr
1242 type(states_elec_t), intent(in) :: st
1243 type(hamiltonian_elec_t), intent(in) :: hm
1244 type(v_ks_t), intent(in) :: ks
1245 type(partner_list_t), intent(in) :: ext_partners
1246 integer, intent(in) :: iter
1247 integer, intent(out) :: ierr
1248
1249 type(gauge_field_t), pointer :: gfield
1250 integer :: err, err2
1251
1252 push_sub(td_dump)
1253
1254 ierr = 0
1255
1256 if (td%restart_dump%skip()) then
1257 pop_sub(td_dump)
1258 return
1259 end if
1260
1261 message(1) = "Debug: Writing td restart."
1262 call messages_info(1, namespace=namespace, debug_only=.true.)
1263
1264 ! first write resume file
1265 call states_elec_dump(td%restart_dump, space, st, gr, hm%kpoints, err, iter=iter)
1266 if (err /= 0) ierr = ierr + 1
1267
1268 call states_elec_dump_rho(td%restart_dump, st, gr, ierr, iter=iter)
1269 if (err /= 0) ierr = ierr + 1
1270
1271 if (hm%lda_u_level /= dft_u_none) then
1272 call lda_u_dump(td%restart_dump, namespace, hm%lda_u, st, gr, ierr)
1273 if (err /= 0) ierr = ierr + 1
1274 end if
1275
1276 call potential_interpolation_dump(td%tr%vks_old, td%restart_dump, gr, st%d%nspin, err2)
1277 if (err2 /= 0) ierr = ierr + 2
1278
1279 call pes_dump(td%pesv, namespace, td%restart_dump, st, gr, err)
1280 if (err /= 0) ierr = ierr + 4
1281
1282 ! Gauge field restart
1283 gfield => list_get_gauge_field(ext_partners)
1284 if(associated(gfield)) then
1285 call gauge_field_dump(td%restart_dump, gfield, ierr)
1286 end if
1287
1288 if (gr%der%boundaries%spiralBC) then
1289 call states_elec_dump_spin(td%restart_dump, st, err)
1290 if (err /= 0) ierr = ierr + 8
1291 end if
1292
1293 if (ks%has_photons) then
1294 call mf_photons_dump(td%restart_dump, ks%pt_mx, gr, td%dt, ks%pt, err)
1295 if (err /= 0) ierr = ierr + 16
1296 end if
1297
1298 if (ks%xc_photon /= 0) then
1299 ! photon-free mean field
1300 call ks%xc_photons%mf_dump(td%restart_dump, err)
1301 if (err /= 0) ierr = ierr + 32
1302 end if
1303
1304 if (allocated(st%frozen_rho)) then
1305 call states_elec_dump_frozen(td%restart_dump, space, st, gr, ierr)
1306 end if
1307 if (err /= 0) ierr = ierr + 64
1308
1309 if (td%ions_dyn%ions_move() .or. td%ions_dyn%cell_relax()) then
1310 call ion_dynamics_dump(td%ions_dyn, td%restart_dump, err)
1311 end if
1312 if (err /= 0) ierr = ierr + 128
1313
1314 message(1) = "Debug: Writing td restart done."
1315 call messages_info(1, namespace=namespace, debug_only=.true.)
1316
1317 pop_sub(td_dump)
1318 end subroutine td_dump
1319
1320 ! ---------------------------------------------------------
1321 subroutine td_load(restart, namespace, space, gr, st, hm, ext_partners, td, ks, ierr)
1322 type(restart_t), intent(in) :: restart
1323 type(namespace_t), intent(in) :: namespace
1324 class(space_t), intent(in) :: space
1325 type(grid_t), intent(in) :: gr
1326 type(states_elec_t), intent(inout) :: st
1327 type(hamiltonian_elec_t), intent(inout) :: hm
1328 type(partner_list_t), intent(in) :: ext_partners
1329 type(td_t), intent(inout) :: td
1330 type(v_ks_t), intent(inout) :: ks
1331 integer, intent(out) :: ierr
1333 integer :: err, err2
1334 type(gauge_field_t), pointer :: gfield
1335 push_sub(td_load)
1336
1337 ierr = 0
1338
1339 if (restart%skip()) then
1340 ierr = -1
1341 pop_sub(td_load)
1342 return
1343 end if
1344
1345 message(1) = "Debug: Reading td restart."
1346 call messages_info(1, namespace=namespace, debug_only=.true.)
1347
1348 ! Read states
1349 call states_elec_load(restart, namespace, space, st, gr, hm%kpoints, err, iter=td%iter, label = ": td")
1350 if (err /= 0) then
1351 ierr = ierr + 1
1352 end if
1353
1354 ! read potential from previous interactions
1355 call potential_interpolation_load(td%tr%vks_old, namespace, restart, gr, st%d%nspin, err2)
1356 if (err2 /= 0) ierr = ierr + 2
1357
1358 if (hm%lda_u_level /= dft_u_none) then
1359 call lda_u_load(restart, hm%lda_u, st, hm%energy%dft_u, ierr)
1360 if (err /= 0) ierr = ierr + 1
1361 end if
1362
1363
1364 ! read PES restart
1365 if (td%pesv%calc_spm .or. td%pesv%calc_mask .or. td%pesv%calc_flux) then
1366 call pes_load(td%pesv, namespace, restart, st, err)
1367 if (err /= 0) ierr = ierr + 4
1368 end if
1369
1370 ! Gauge field restart
1371 gfield => list_get_gauge_field(ext_partners)
1372 if (associated(gfield)) then
1373 call gauge_field_load(restart, gfield, err)
1374 if (err /= 0) then
1375 ierr = ierr + 8
1376 else
1377 call hm%update(gr, namespace, space, ext_partners, time = td%dt*td%iter)
1378 end if
1379 end if
1380
1381 ! add photon restart
1382 if (ks%has_photons) then
1383 call mf_photons_load(restart, ks%pt_mx, gr, err)
1384 end if
1385 if (err /= 0) ierr = ierr + 16
1386
1387 if (ks%xc_photon /= 0) then
1388 call ks%xc_photons%mf_load(restart, space, err)
1389 end if
1390 if (err /= 0) ierr = ierr + 32
1391
1392 if (gr%der%boundaries%spiralBC) then
1393 call states_elec_load_spin(restart, st, err)
1394 !To ensure back compatibility, if the file is not present, we use the
1395 !current states to get the spins
1396 if (err /= 0) call states_elec_fermi(st, namespace, gr)
1397 end if
1398
1399 if (td%ions_dyn%is_active()) then
1400 call ion_dynamics_load(td%ions_dyn, restart, err)
1401 end if
1402 if (err /= 0) ierr = ierr + 64
1403
1404 message(1) = "Debug: Reading td restart done."
1405 call messages_info(1, namespace=namespace, debug_only=.true.)
1406
1407 pop_sub(td_load)
1408 end subroutine td_load
1409 ! ---------------------------------------------------------
1410 subroutine td_load_frozen(namespace, restart, space, mesh, st, hm, ierr)
1411 type(namespace_t), intent(in) :: namespace
1412 type(restart_t), intent(in) :: restart
1413 class(space_t), intent(in) :: space
1414 class(mesh_t), intent(in) :: mesh
1415 type(states_elec_t), intent(inout) :: st
1416 type(hamiltonian_elec_t), intent(inout) :: hm
1417 integer, intent(out) :: ierr
1418
1419 push_sub(td_load_frozen)
1420
1421 ierr = 0
1422
1423 if (restart%skip()) then
1424 ierr = -1
1425 pop_sub(td_load_frozen)
1426 return
1427 end if
1428
1429 message(1) = "Debug: Reading td frozen restart."
1430 call messages_info(1, namespace=namespace, debug_only=.true.)
1431
1432 safe_allocate(st%frozen_rho(1:mesh%np, 1:st%d%nspin))
1433 if (family_is_mgga(hm%xc%family)) then
1434 safe_allocate(st%frozen_tau(1:mesh%np, 1:st%d%nspin))
1435 safe_allocate(st%frozen_gdens(1:mesh%np, 1:space%dim, 1:st%d%nspin))
1436 safe_allocate(st%frozen_ldens(1:mesh%np, 1:st%d%nspin))
1437 end if
1438
1439 call states_elec_load_frozen(restart, space, st, mesh, ierr)
1440
1441 message(1) = "Debug: Reading td frozen restart done."
1442 call messages_info(1, namespace=namespace, debug_only=.true.)
1443
1444 pop_sub(td_load_frozen)
1445 end subroutine td_load_frozen
1446
1447 ! ---------------------------------------------------------
1448 logical function td_get_from_scratch(td)
1449 type(td_t), intent(in) :: td
1450
1451 push_sub(td_get_from_scratch)
1452
1453 td_get_from_scratch = td%from_scratch
1454
1455 pop_sub(td_get_from_scratch)
1456 end function td_get_from_scratch
1457
1458 ! ---------------------------------------------------------
1459 subroutine td_set_from_scratch(td, from_scratch)
1460 type(td_t), intent(inout) :: td
1461 logical, intent(in) :: from_scratch
1462
1463 push_sub(td_set_from_scratch)
1464
1465 td%from_scratch = from_scratch
1466
1467 pop_sub(td_set_from_scratch)
1468 end subroutine td_set_from_scratch
1469end module td_oct_m
1470
1471!! Local Variables:
1472!! mode: f90
1473!! coding: utf-8
1474!! End:
Define which routines can be seen from the outside.
Definition: loct.F90:149
Prints out to iunit a message in the form: ["InputVariable" = value] where "InputVariable" is given b...
Definition: messages.F90:182
integer, parameter, public mask_absorbing
Module implementing boundary conditions in Octopus.
Definition: boundaries.F90:124
This module handles the calculation mode.
This module implements a calculator for the density and defines related functions.
Definition: density.F90:122
subroutine, public states_elec_freeze_adjust_qtot(st)
Definition: density.F90:814
subroutine, public states_elec_freeze_redistribute_states(st, namespace, mesh, mc, nn)
Definition: density.F90:786
subroutine, public states_elec_freeze_orbitals(st, namespace, space, gr, mc, kpoints, n, family_is_mgga)
Calculate partial density for frozen orbitals.
Definition: density.F90:655
subroutine, public density_calc(st, gr, density, istin)
Computes the density from the orbitals in st.
Definition: density.F90:612
subroutine, public dm_propagation_init_run(adiabatic_st, namespace, space, gr, st, hm, mc)
Initialise the adiabatic states prior to running TB propagation.
subroutine, public dm_propagation_run(dmp, namespace, space, gr, ions, st, mc, hm, ks, iter, dt, ext_partners, update_energy)
Apply dissipation to a TD run via the Linblad formalism.
A set of subroutines for performing the parts of a ground state calculation with an electrons system....
subroutine, public electrons_ground_state_run(namespace, mc, gr, ions, ext_partners, st, ks, hm, outp, space, fromScratch)
Run a ground state calculation for a system of electrons.
subroutine, public energy_calc_total(namespace, space, hm, gr, st, ext_partners, iunit, full)
This subroutine calculates the total energy of the system. Basically, it adds up the KS eigenvalues,...
integer, parameter, public spin_orbit
Definition: epot.F90:168
logical function, public list_has_gauge_field(partners)
type(gauge_field_t) function, pointer, public list_get_gauge_field(partners)
type(lasers_t) function, pointer, public list_get_lasers(partners)
subroutine, public forces_calculate(gr, namespace, ions, hm, ext_partners, st, ks, vhxc_old, t, dt)
Definition: forces.F90:340
subroutine, public gauge_field_load(restart, gfield, ierr)
subroutine, public gauge_field_get_force(this, gr, spin_channels, current, lrc)
subroutine, public gauge_field_dump(restart, gfield, ierr)
logical pure function, public gauge_field_is_propagated(this)
logical pure function, public gauge_field_is_used(this)
subroutine, public gauge_field_init_vec_pot(this, qtot)
real(real64), parameter, public m_zero
Definition: global.F90:191
real(real64), parameter, public m_epsilon
Definition: global.F90:207
real(real64), parameter, public m_one
Definition: global.F90:192
This module implements the underlying real-space grid.
Definition: grid.F90:119
subroutine, public zvmask(mesh, hm, st)
subroutine, public hamiltonian_elec_epot_generate(this, namespace, space, gr, ions, ext_partners, st, time)
subroutine, public dvmask(mesh, hm, st)
This module defines classes and functions for interaction partners.
Definition: io.F90:116
subroutine, public io_close(iunit, grp)
Definition: io.F90:466
subroutine, public io_skip_header(iunit)
Definition: io.F90:645
character(len=max_path_len) function, public io_workpath(path, namespace)
construct path name from given name and namespace
Definition: io.F90:317
subroutine, public io_debug_on_the_fly(namespace)
check if debug mode should be enabled or disabled on the fly
Definition: io.F90:534
integer function, public io_open(file, namespace, action, status, form, position, die, recl, grp)
Definition: io.F90:401
subroutine, public ion_dynamics_dump(this, restart, ierr)
subroutine, public ion_dynamics_propagate(this, ions, time, dt, namespace)
Interface for the ion/cell dynamics.
subroutine, public ion_dynamics_load(this, restart, ierr)
subroutine, public ion_dynamics_init(this, namespace, ions, symmetrize, symm)
subroutine, public ion_dynamics_end(this)
logical pure function, public ion_dynamics_drive_ions(this)
Is the ion dynamics activated or not.
subroutine, public kick_apply(space, mesh, st, ions_dyn, ions, kick, psolver, kpoints, pcm)
Applies the delta-function electric field where k = kick%delta_strength.
Definition: kick.F90:1137
subroutine, public lasers_set_nondipole_parameters(this, ndfield, nd_integration_time)
Set parameters for nondipole SFA calculation.
Definition: lasers.F90:754
logical function, public lasers_with_nondipole_field(lasers)
Check if a nondipole SFA correction should be computed for the given laser.
Definition: lasers.F90:741
subroutine, public laser_write_info(lasers, namespace, dt, max_iter, iunit)
Definition: lasers.F90:932
subroutine, public lda_u_dump(restart, namespace, this, st, mesh, ierr)
Definition: lda_u_io.F90:648
subroutine, public lda_u_write_u(this, iunit, namespace)
Definition: lda_u_io.F90:532
subroutine, public lda_u_load(restart, this, st, dftu_energy, ierr, occ_only, u_only)
Definition: lda_u_io.F90:729
subroutine, public lda_u_write_v(this, iunit, namespace)
Definition: lda_u_io.F90:581
integer, parameter, public dft_u_none
Definition: lda_u.F90:203
subroutine, public lda_u_init(this, namespace, space, level, gr, ions, st, mc, kpoints)
Definition: lda_u.F90:284
subroutine, public lda_u_freeze_occ(this)
Definition: lda_u.F90:887
subroutine, public lda_u_freeze_u(this)
Definition: lda_u.F90:894
subroutine, public lda_u_end(this)
Definition: lda_u.F90:655
subroutine, public lda_u_update_occ_matrices(this, namespace, mesh, st, phase, energy)
Definition: lda_u.F90:798
integer, parameter, public dft_u_acbn0
Definition: lda_u.F90:203
This module implements fully polymorphic linked lists, and some specializations thereof.
This module defines the meshes, which are used in Octopus.
Definition: mesh.F90:120
subroutine, public messages_print_with_emphasis(msg, iunit, namespace)
Definition: messages.F90:898
subroutine, public messages_not_implemented(feature, namespace)
Definition: messages.F90:1091
character(len=512), private msg
Definition: messages.F90:167
subroutine, public messages_warning(no_lines, all_nodes, namespace)
Definition: messages.F90:525
subroutine, public messages_obsolete_variable(namespace, name, rep)
Definition: messages.F90:1023
character(len=256), dimension(max_lines), public message
to be output by fatal, warning
Definition: messages.F90:162
subroutine, public messages_fatal(no_lines, only_root_writes, namespace)
Definition: messages.F90:410
subroutine, public messages_input_error(namespace, var, details, row, column)
Definition: messages.F90:691
subroutine, public messages_experimental(name, namespace)
Definition: messages.F90:1063
subroutine, public messages_info(no_lines, iunit, debug_only, stress, all_nodes, namespace)
Definition: messages.F90:594
type(mpi_grp_t), public mpi_world
Definition: mpi.F90:272
This module handles the communicators for the various parallelization strategies.
Definition: multicomm.F90:147
type(namespace_t), public global_namespace
Definition: namespace.F90:134
this module contains the low-level part of the output system
Definition: output_low.F90:117
this module contains the output system
Definition: output.F90:117
logical function, public parse_is_defined(namespace, name)
Definition: parser.F90:455
subroutine, public pes_calc(pes, namespace, space, mesh, st, dt, iter, der, kpoints, ext_partners, stopping)
Definition: pes.F90:271
subroutine, public pes_init(pes, namespace, space, mesh, box, st, save_iter, kpoints, abs_boundaries, ext_partners, max_iter, dt)
Definition: pes.F90:177
subroutine, public pes_output(pes, namespace, space, gr, st, iter, outp, dt, ions)
Definition: pes.F90:297
subroutine, public pes_init_write(pes, mesh, st, namespace)
Definition: pes.F90:401
subroutine, public pes_end(pes)
Definition: pes.F90:257
subroutine, public pes_load(pes, namespace, restart, st, ierr)
Definition: pes.F90:361
subroutine, public pes_dump(pes, namespace, restart, st, mesh, ierr)
Definition: pes.F90:321
subroutine, public mf_photons_load(restart, this, gr, ierr)
subroutine, public mf_photons_dump(restart, this, gr, dt, pt_mode, ierr)
subroutine, public potential_interpolation_load(potential_interpolation, namespace, restart, mesh, nspin, err2)
subroutine, public potential_interpolation_dump(potential_interpolation, restart, mesh, nspin, err2)
subroutine, public profiling_out(label)
Increment out counter and sum up difference between entry and exit time.
Definition: profiling.F90:625
subroutine, public profiling_in(label, exclude)
Increment in counter and save entry time.
Definition: profiling.F90:554
subroutine, public propagator_elec_dt_bo(scf, namespace, space, gr, ks, st, hm, ions, ext_partners, mc, iter, dt, ions_dyn, scsteps)
subroutine, public propagator_elec_init(gr, namespace, st, tr, ks_pot, have_fields, family_is_mgga_with_exc, relax_cell)
subroutine, public propagator_elec_dt(ks, namespace, space, hm, gr, st, tr, time, dt, nt, ions_dyn, ions, ext_partners, mc, outp, write_handler, scsteps, update_energy, qcchi)
Propagates st from time - dt to t. If dt<0, it propagates backwards from t+|dt| to t.
subroutine, public propagator_elec_end(tr)
This module implements the basic propagator framework.
Definition: propagator.F90:119
logical function, public clean_stop(comm)
returns true if a file named stop exists
Definition: restart.F90:335
integer, parameter, public restart_gs
Definition: restart.F90:156
integer, parameter, public restart_type_dump
Definition: restart.F90:183
integer, parameter, public restart_td
Definition: restart.F90:156
integer, parameter, public restart_type_load
Definition: restart.F90:183
subroutine, public scf_init(scf, namespace, gr, ions, st, mc, hm, space)
Definition: scf.F90:256
subroutine, public scf_end(scf)
Definition: scf.F90:526
subroutine, public scissor_init(this, namespace, space, st, mesh, kpoints, phase, gap, mc)
Definition: scissor.F90:164
pure logical function, public states_are_real(st)
subroutine, public states_elec_fermi(st, namespace, mesh, compute_spin)
calculate the Fermi level for the states in this object
subroutine, public states_elec_end(st)
finalize the states_elec_t object
subroutine, public states_elec_deallocate_wfns(st)
Deallocates the KS wavefunctions defined within a states_elec_t structure.
subroutine, public states_elec_allocate_wfns(st, mesh, wfs_type, skip, packed)
Allocates the KS wavefunctions defined within a states_elec_t structure.
This module handles reading and writing restart information for the states_elec_t.
subroutine, public states_elec_read_user_def_orbitals(mesh, namespace, space, st)
the routine reads formulas for user-defined wavefunctions from the input file and fills the respectiv...
subroutine, public states_elec_load_frozen(restart, space, st, mesh, ierr)
subroutine, public states_elec_transform(st, namespace, space, restart, mesh, kpoints, prefix)
subroutine, public states_elec_dump(restart, space, st, mesh, kpoints, ierr, iter, lr, verbose)
subroutine, public states_elec_load(restart, namespace, space, st, mesh, kpoints, ierr, iter, lr, lowest_missing, label, verbose, skip)
returns in ierr: <0 => Fatal error, or nothing read =0 => read all wavefunctions >0 => could only rea...
subroutine, public states_elec_load_spin(restart, st, ierr)
returns in ierr: <0 => Fatal error, or nothing read =0 => read all wavefunctions >0 => could only rea...
subroutine, public states_elec_dump_frozen(restart, space, st, mesh, ierr)
subroutine, public states_elec_dump_rho(restart, st, mesh, ierr, iter)
subroutine, public states_elec_dump_spin(restart, st, ierr)
This module implements the calculation of the stress tensor.
Definition: stress.F90:120
subroutine, public stress_calculate(namespace, gr, hm, st, ions, ks, ext_partners)
This computes the total stress on the lattice.
Definition: stress.F90:188
Definition: td.F90:116
subroutine, public td_check_point(td, namespace, mc, gr, ions, st, ks, hm, ext_partners, outp, space, iter, scsteps, etime, stopping, from_scratch)
Definition: td.F90:793
subroutine, public td_end(td)
Definition: td.F90:632
subroutine, public td_load_restart_from_gs(td, namespace, space, mc, gr, ext_partners, st, ks, hm)
Definition: td.F90:1202
subroutine td_print_message(td, namespace, ions, hm, iter, scsteps, etime)
Definition: td.F90:857
subroutine, public td_end_run(td, st, hm)
Definition: td.F90:651
subroutine, public td_init(td, namespace, space, gr, ions, st, ks, hm, ext_partners, outp)
Definition: td.F90:240
subroutine, public td_allocate_wavefunctions(td, namespace, mc, gr, ions, st, hm, space)
Definition: td.F90:568
logical function, public td_get_from_scratch(td)
Definition: td.F90:1544
subroutine td_load_frozen(namespace, restart, space, mesh, st, hm, ierr)
Definition: td.F90:1506
subroutine td_print_header(namespace)
Definition: td.F90:778
integer, parameter, public bo
Definition: td.F90:203
subroutine td_load(restart, namespace, space, gr, st, hm, ext_partners, td, ks, ierr)
Definition: td.F90:1417
subroutine, public td_set_from_scratch(td, from_scratch)
Definition: td.F90:1555
subroutine, public td_dump(td, namespace, space, gr, st, hm, ks, ext_partners, iter, ierr)
Definition: td.F90:1333
subroutine td_read_coordinates(td, namespace, ions)
reads the pos and vel from coordinates file
Definition: td.F90:1290
subroutine, public td_run(td, namespace, mc, gr, ions, st, ks, hm, ext_partners, outp, space, from_scratch)
Definition: td.F90:671
subroutine td_run_zero_iter(td, namespace, space, gr, ions, st, ks, hm, ext_partners, outp, mc)
Definition: td.F90:1244
subroutine, public td_init_gaugefield(td, namespace, gr, st, ks, hm, ext_partners, space)
Definition: td.F90:601
subroutine td_update_elapsed_time(etime)
Definition: td.F90:878
subroutine, public td_init_with_wavefunctions(td, namespace, space, mc, gr, ions, ext_partners, st, ks, hm, outp, from_scratch)
Definition: td.F90:889
subroutine, public td_load_restart_from_td(td, namespace, space, mc, gr, ext_partners, st, ks, hm, from_scratch)
Definition: td.F90:1166
subroutine td_init_ions_and_forces(td, namespace, space, gr, ions, ext_partners, st, ks, hm, outp)
Definition: td.F90:1120
subroutine, public td_init_run(td, namespace, mc, gr, ions, st, ks, hm, ext_partners, outp, space, from_scratch)
Definition: td.F90:512
subroutine, public td_write_output(namespace, space, gr, st, hm, ks, outp, ions, ext_partners, iter, dt)
Definition: td_write.F90:1239
subroutine, public td_write_data(writ)
Definition: td_write.F90:1205
subroutine, public td_write_kick(outp, namespace, space, mesh, kick, ions, iter)
Definition: td_write.F90:337
subroutine, public td_write_end(writ)
Definition: td_write.F90:989
subroutine, public td_write_iter(writ, namespace, space, outp, gr, st, hm, ions, ext_partners, kick, ks, dt, iter, mc, recalculate_gs)
Definition: td_write.F90:1036
subroutine, public td_write_init(writ, namespace, space, outp, gr, st, hm, ions, ext_partners, ks, ions_move, with_gauge_field, kick, iter, max_iter, dt, mc)
Initialize files to write when prograting in time.
Definition: td_write.F90:375
integer, parameter, public out_separate_forces
Definition: td_write.F90:203
type(type_t), public type_cmplx
Definition: types.F90:136
brief This module defines the class unit_t which is used by the unit_systems_oct_m module.
Definition: unit.F90:134
This module defines the unit system, used for input and output.
type(unit_t), public unit_femtosecond
Time in femtoseconds.
type(unit_system_t), public units_out
type(unit_system_t), public units_inp
the units systems for reading and writing
subroutine, public v_ks_freeze_hxc(ks)
Definition: v_ks.F90:1483
subroutine, public v_ks_calculate_current(this, calc_cur)
Definition: v_ks.F90:1494
subroutine, public v_ks_calc(ks, namespace, space, hm, st, ions, ext_partners, calc_eigenval, time, calc_energy, calc_current, force_semilocal)
Definition: v_ks.F90:749
This module provices a simple timer class which can be used to trigger the writing of a restart file ...
Definition: walltimer.F90:123
logical function, public walltimer_alarm(comm, print)
indicate whether time is up
Definition: walltimer.F90:333
real(real64) function, public walltimer_get_start_time()
Return the walltimer start time.
Definition: walltimer.F90:407
Definition: xc.F90:116
pure logical function, public family_is_mgga(family, only_collinear)
Is the xc function part of the mGGA family.
Definition: xc.F90:580
logical pure function, public family_is_mgga_with_exc(xcs)
Is the xc function part of the mGGA family with an energy functional.
Definition: xc.F90:593
Extension of space that contains the knowledge of the spin dimension.
Description of the grid, containing information on derivatives, stencil, and symmetries.
Definition: grid.F90:171
Describes mesh distribution to nodes.
Definition: mesh.F90:187
Stores all communicators and groups.
Definition: multicomm.F90:208
output handler class
Definition: output_low.F90:166
The states_elec_t class contains all electronic wave functions.
int true(void)