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
594 push_sub(td_run)
595
596 etime = loct_clock()
597 ! This is the time-propagation loop. It starts at t=0 and finishes
598 ! at td%max_iter*dt. The index i runs from 1 to td%max_iter, and
599 ! step "iter" means propagation from (iter-1)*dt to iter*dt.
600 propagation: do iter = td%iter, td%max_iter
601
602 stopping = clean_stop(mc%master_comm) .or. walltimer_alarm(mc%master_comm)
603
604 call profiling_in("TIME_STEP")
605
606 if (iter > 1) then
607 if (((iter-1)*td%dt <= hm%kick%time) .and. (iter*td%dt > hm%kick%time)) then
608 if (.not. hm%pcm%localf) then
609 call kick_apply(space, gr, st, td%ions_dyn, ions, hm%kick, hm%psolver, hm%kpoints)
610 else
611 call kick_apply(space, gr, st, td%ions_dyn, ions, hm%kick, hm%psolver, hm%kpoints, pcm = hm%pcm)
612 end if
613 call td_write_kick(outp, namespace, space, gr, hm%kick, ions, iter)
614 !We activate the sprial BC only after the kick,
615 !to be sure that the first iteration corresponds to the ground state
616 if (gr%der%boundaries%spiralBC) gr%der%boundaries%spiral = .true.
617 end if
618 end if
619
620 ! time iterate the system, one time step.
621 select case (td%dynamics)
622 case (ehrenfest)
623 call propagator_elec_dt(ks, namespace, space, hm, gr, st, td%tr, iter*td%dt, td%dt, iter, td%ions_dyn, &
624 ions, ext_partners, mc, outp, td%write_handler, scsteps = scsteps, &
625 update_energy = (mod(iter, td%energy_update_iter) == 0) .or. (iter == td%max_iter))
626 case (bo)
627 call propagator_elec_dt_bo(td%scf, namespace, space, gr, ks, st, hm, ions, ext_partners, mc, iter, &
628 td%dt, td%ions_dyn, scsteps)
629 end select
630
631 if (td%dmp%calculation_mode /= option__tddmpropagation__no_propagation) then
632 call dm_propagation_run(td%dmp, namespace, space, gr, ions, st, mc, hm, ks, iter, td%dt, ext_partners, &
633 update_energy = (mod(iter, td%energy_update_iter) == 0) .or. (iter == td%max_iter))
634 end if
635
636 !Apply mask absorbing boundaries
637 if (hm%abs_boundaries%abtype == mask_absorbing) then
638 if (states_are_real(st)) then
639 call dvmask(gr, hm, st)
640 else
641 call zvmask(gr, hm, st)
642 end if
643 end if
644
645 !Photoelectron stuff
646 if (td%pesv%calc_spm .or. td%pesv%calc_mask .or. td%pesv%calc_flux) then
647 call pes_calc(td%pesv, namespace, space, gr, st, td%dt, iter, gr%der, hm%kpoints, ext_partners, stopping)
648 end if
649
650 call td_write_iter(td%write_handler, namespace, space, outp, gr, st, hm, ions, ext_partners, &
651 hm%kick, ks, td%dt, iter, mc, td%recalculate_gs)
652
653 ! write down data
654 call td_check_point(td, namespace, mc, gr, ions, st, ks, hm, ext_partners, outp, space, &
655 iter, scsteps, etime, stopping, from_scratch)
656
657 ! check if debug mode should be enabled or disabled on the fly
658 call io_debug_on_the_fly(namespace)
659
660 call profiling_out("TIME_STEP")
661 if (stopping) exit
662
663 end do propagation
664
665 ! Output propagation speed in fs/day
666 wall_time = loct_clock() - walltimer_get_start_time()
667 simulation_time = td%dt * (iter - td%iter + 1)
668 speed_fs_per_day = units_from_atomic(unit_femtosecond, simulation_time) / (wall_time / 86400.0_real64)
669 write(message(1), '(a,f10.3,a)') 'Propagation speed: ', speed_fs_per_day, ' fs/day'
670 write(message(2), '(a)') ''
671 call messages_info(2, namespace=namespace)
672
673 pop_sub(td_run)
674 end subroutine td_run
675
676 subroutine td_print_header(namespace)
677 type(namespace_t), intent(in) :: namespace
678
679 push_sub(td_print_header)
680
681 write(message(1), '(a7,1x,a14,a14,a10,a17)') 'Iter ', 'Time ', 'Energy ', 'SC Steps', 'Elapsed Time '
682
683 call messages_info(1, namespace=namespace)
684 call messages_print_with_emphasis(namespace=namespace)
685
686 pop_sub(td_print_header)
687 end subroutine td_print_header
688
689 ! ---------------------------------------------------------
690 subroutine td_check_point(td, namespace, mc, gr, ions, st, ks, hm, ext_partners, outp, space, &
691 iter, scsteps, etime, stopping, from_scratch)
692 type(td_t), intent(inout) :: td
693 type(namespace_t), intent(in) :: namespace
694 type(multicomm_t), intent(in) :: mc
695 type(grid_t), intent(inout) :: gr
696 type(ions_t), intent(inout) :: ions
697 type(states_elec_t), intent(inout) :: st
698 type(v_ks_t), intent(inout) :: ks
699 type(hamiltonian_elec_t), intent(inout) :: hm
700 type(partner_list_t), intent(in) :: ext_partners
701 type(output_t), intent(in) :: outp
702 type(electron_space_t), intent(in) :: space
703 integer, intent(in) :: iter
704 integer, intent(in) :: scsteps
705 real(real64), intent(inout) :: etime
706 logical, intent(in) :: stopping
707 logical, intent(inout) :: from_scratch
708
709 integer :: ierr
710
711 push_sub(td_check_point)
712
713 call td_print_message(td, namespace, ions, hm, iter, scsteps, etime)
714
715 if (outp%anything_now(iter)) then ! output
716 call td_write_output(namespace, space, gr, st, hm, ks, outp, ions, ext_partners, iter, td%dt)
717 end if
718
719 if (mod(iter, outp%restart_write_interval) == 0 .or. iter == td%max_iter .or. stopping) then ! restart
720 !if (iter == td%max_iter) outp%iter = ii - 1
721 call td_write_data(td%write_handler)
722 call td_dump(td, namespace, space, gr, st, hm, ks, ext_partners, iter, ierr)
723 if (ierr /= 0) then
724 message(1) = "Unable to write time-dependent restart information."
725 call messages_warning(1, namespace=namespace)
726 end if
727
728 call pes_output(td%pesv, namespace, space, gr, st, iter, outp, td%dt, ions)
729
730 if ((td%ions_dyn%is_active()) .and. td%recalculate_gs) then
731 call messages_print_with_emphasis(msg='Recalculating the ground state.', namespace=namespace)
732 from_scratch = .false.
734 call electrons_ground_state_run(namespace, mc, gr, ions, ext_partners, st, ks, hm, outp, space, from_scratch)
735 call states_elec_allocate_wfns(st, gr, packed=.true.)
736 call td_load(td%restart_load, namespace, space, gr, st, hm, ext_partners, td, ks, ierr)
737 if (ierr /= 0) then
738 message(1) = "Unable to load TD states."
739 call messages_fatal(1, namespace=namespace)
740 end if
741 call density_calc(st, gr, st%rho)
742 call v_ks_calc(ks, namespace, space, hm, st, ions, ext_partners, &
743 calc_eigenval=.true., time = iter*td%dt, calc_energy=.true.)
744 call forces_calculate(gr, namespace, ions, hm, ext_partners, st, ks, t = iter*td%dt, dt = td%dt)
745 assert(.not. td%ions_dyn%cell_relax())
746 call messages_print_with_emphasis(msg="Time-dependent simulation proceeds", namespace=namespace)
747 call td_print_header(namespace)
748 end if
749 end if
750
751 pop_sub(td_check_point)
752 end subroutine td_check_point
753
754 ! ---------------------------------------------------------
755 subroutine td_print_message(td, namespace, ions, hm, iter, scsteps, etime)
756 type(td_t), intent(inout) :: td
757 type(namespace_t), intent(in) :: namespace
758 type(ions_t), intent(inout) :: ions
759 type(hamiltonian_elec_t), intent(inout) :: hm
760 integer, intent(in) :: iter
761 integer, intent(in) :: scsteps
762 real(real64), intent(inout) :: etime
763
764 push_sub(td_print_message)
765
766 write(message(1), '(i7,1x,2f14.6,i10,f14.3)') iter, units_from_atomic(units_out%time, iter*td%dt), &
767 units_from_atomic(units_out%energy, hm%energy%total + ions%kinetic_energy), scsteps, &
768 loct_clock() - etime
769 call messages_info(1, namespace=namespace)
770 call td_update_elapsed_time(etime)
772 pop_sub(td_print_message)
773 end subroutine td_print_message
774
775 ! ---------------------------------------------------------
776 subroutine td_update_elapsed_time(etime)
777 real(real64), intent(inout) :: etime
778
779 push_sub(td_update_elapsed_time)
780
781 etime = loct_clock()
782
784 end subroutine td_update_elapsed_time
786 ! ---------------------------------------------------------
787 subroutine td_init_with_wavefunctions(td, namespace, space, mc, gr, ions, ext_partners, st, ks, hm, outp, from_scratch)
788 type(td_t), intent(inout) :: td
789 type(namespace_t), intent(in) :: namespace
790 type(electron_space_t), intent(in) :: space
791 type(multicomm_t), intent(in) :: mc
792 type(grid_t), intent(inout) :: gr
793 type(ions_t), intent(inout) :: ions
794 type(partner_list_t), intent(in) :: ext_partners
795 type(states_elec_t), target, intent(inout) :: st
796 type(v_ks_t), intent(inout) :: ks
797 type(hamiltonian_elec_t), intent(inout) :: hm
798 type(output_t), intent(inout) :: outp
799 logical, intent(in) :: from_scratch
800
801 integer :: ierr
802 real(real64) :: x
803 real(real64) :: ndinitial(space%dim)
804 logical :: freeze_hxc, freeze_occ, freeze_u
805 type(restart_t) :: restart, restart_frozen
806 type(gauge_field_t), pointer :: gfield
807 type(lasers_t), pointer :: lasers
809
810 !We activate the sprial BC only after the kick,
811 !to be sure that the first iteration corresponds to the ground state
812 if (gr%der%boundaries%spiralBC) then
813 if ((td%iter-1)*td%dt > hm%kick%time) then
814 gr%der%boundaries%spiral = .true.
815 end if
816 hm%vnl%spin => st%spin
817 hm%phase%spin => st%spin
818 !We fill st%spin. In case of restart, we read it in td_load
819 if (from_scratch) call states_elec_fermi(st, namespace, gr)
820 end if
821
822 if (from_scratch) then
823 ! Initialize the occupation matrices and U for DFT+U
824 ! This must be called before parsing TDFreezeOccupations and TDFreezeU
825 ! in order that the code does properly the initialization.
826 call lda_u_update_occ_matrices(hm%lda_u, namespace, gr, st, hm%phase, hm%energy)
827 end if
828
829 if (td%freeze_orbitals > 0) then
830 if (from_scratch) then
831 ! In this case, we first freeze the orbitals, then calculate the Hxc potential.
832 call states_elec_freeze_orbitals(st, namespace, space, gr, mc, hm%kpoints, &
833 td%freeze_orbitals, family_is_mgga(ks%xc_family))
834 else
835 call restart%init(namespace, restart_td, restart_type_load, mc, ierr, mesh=gr)
836 if (ierr == 0) then
837 call td_load_frozen(namespace, restart, space, gr, st, hm, ierr)
838 end if
839 if (ierr /= 0) then
840 td%iter = 0
841 message(1) = "Unable to read frozen restart information."
842 call messages_fatal(1, namespace=namespace)
843 end if
844 call restart%end()
845 end if
846 write(message(1),'(a,i4,a,i4,a)') 'Info: The lowest', td%freeze_orbitals, &
847 ' orbitals have been frozen.', st%nst, ' will be propagated.'
848 call messages_info(1, namespace=namespace)
850 call density_calc(st, gr, st%rho)
851 call v_ks_calc(ks, namespace, space, hm, st, ions, ext_partners, calc_eigenval=.true., time = td%iter*td%dt)
852 else if (td%freeze_orbitals < 0) then
853 ! This means SAE approximation. We calculate the Hxc first, then freeze all
854 ! orbitals minus one.
855 write(message(1),'(a)') 'Info: The single-active-electron approximation will be used.'
856 call messages_info(1, namespace=namespace)
857 call v_ks_calc(ks, namespace, space, hm, st, ions, ext_partners, calc_eigenval=.true., time = td%iter*td%dt)
858 if (from_scratch) then
859 call states_elec_freeze_orbitals(st, namespace, space, gr, mc, hm%kpoints, st%nst-1, family_is_mgga(ks%xc_family))
860 else
861 call messages_not_implemented("TDFreezeOrbials < 0 with FromScratch=no", namespace=namespace)
862 end if
864 call v_ks_freeze_hxc(ks)
865 call density_calc(st, gr, st%rho)
866 else
867 ! Normal run.
868 call density_calc(st, gr, st%rho)
869 call v_ks_calc(ks, namespace, space, hm, st, ions, ext_partners, calc_eigenval=.true., time = td%iter*td%dt)
870 end if
872 !%Variable TDFreezeHXC
873 !%Type logical
874 !%Default no
875 !%Section Time-Dependent
876 !%Description
877 !% The electrons are evolved as independent particles feeling the Hartree and
878 !% exchange-correlation potentials from the ground-state electronic configuration.
879 !%End
880 call parse_variable(namespace, 'TDFreezeHXC', .false., freeze_hxc)
881 if (freeze_hxc) then
882 write(message(1),'(a)') 'Info: Freezing Hartree and exchange-correlation potentials.'
883 call messages_info(1, namespace=namespace)
884
885 if (.not. from_scratch) then
886
887 call restart_frozen%init(namespace, restart_gs, restart_type_load, mc, ierr, mesh=gr, exact=.true.)
888 call states_elec_load(restart_frozen, namespace, space, st, gr, hm%kpoints, ierr, label = ": gs")
889 call states_elec_transform(st, namespace, space, restart_frozen, gr, hm%kpoints)
890 call restart_frozen%end()
891
892 call density_calc(st, gr, st%rho)
893 call v_ks_calc(ks, namespace, space, hm, st, ions, ext_partners, calc_eigenval=.true., time = td%iter*td%dt)
894
895 call restart_frozen%init(namespace, restart_td, restart_type_load, mc, ierr, mesh=gr)
896 call states_elec_load(restart_frozen, namespace, space, st, gr, hm%kpoints, ierr, iter=td%iter, label = ": td")
897 call restart_frozen%end()
898 call hm%ks_pot%run_zero_iter(td%tr%vks_old)
899
900 end if
901
902 call v_ks_freeze_hxc(ks)
903
904 end if
905
906 x = minval(st%eigenval(st%st_start, :))
907 if (st%parallel_in_states) then
908 call st%mpi_grp%bcast(x, 1, mpi_double_precision, 0)
909 end if
910 call hm%update_span(gr%spacing(1:space%dim), x, namespace)
911 ! initialize Fermi energy
912 call states_elec_fermi(st, namespace, gr, compute_spin = .not. gr%der%boundaries%spiralBC)
913 call energy_calc_total(namespace, space, hm, gr, st, ext_partners)
914
915 !%Variable TDFreezeDFTUOccupations
916 !%Type logical
917 !%Default no
918 !%Section Time-Dependent
919 !%Description
920 !% The occupation matrices than enters in the DFT+U potential
921 !% are not evolved during the time evolution.
922 !%End
923 call parse_variable(namespace, 'TDFreezeDFTUOccupations', .false., freeze_occ)
924 if (freeze_occ) then
925 write(message(1),'(a)') 'Info: Freezing DFT+U occupation matrices that enters in the DFT+U potential.'
926 call messages_info(1, namespace=namespace)
927 call lda_u_freeze_occ(hm%lda_u)
928
929 !In this case we should reload GS wavefunctions
930 if (hm%lda_u_level /= dft_u_none .and. .not. from_scratch) then
931 call restart_frozen%init(namespace, restart_gs, restart_type_load, mc, ierr, mesh=gr)
932 call lda_u_load(restart_frozen, hm%lda_u, st, hm%energy%dft_u, ierr, occ_only = .true.)
933 call restart_frozen%end()
934 end if
935 end if
936
937 !%Variable TDFreezeU
938 !%Type logical
939 !%Default no
940 !%Section Time-Dependent
941 !%Description
942 !% The effective U of DFT+U is not evolved during the time evolution.
943 !%End
944 call parse_variable(namespace, 'TDFreezeU', .false., freeze_u)
945 if (freeze_u) then
946 write(message(1),'(a)') 'Info: Freezing the effective U of DFT+U.'
947 call messages_info(1, namespace=namespace)
948 call lda_u_freeze_u(hm%lda_u)
949
950 !In this case we should reload GS wavefunctions
951 if (hm%lda_u_level == dft_u_acbn0 .and. .not. from_scratch) then
952 call restart_frozen%init(namespace, restart_gs, restart_type_load, mc, ierr, mesh=gr)
953 call lda_u_load(restart_frozen, hm%lda_u, st, hm%energy%dft_u, ierr, u_only = .true.)
954 call restart_frozen%end()
955 write(message(1),'(a)') 'Loaded GS effective U of DFT+U'
956 call messages_info(1, namespace=namespace)
957 call lda_u_write_u(hm%lda_u, namespace=namespace)
958 call lda_u_write_v(hm%lda_u, namespace=namespace)
959 end if
960 end if
961
962 ! This needs to be called before the calculation of the forces,
963 ! as we need to test of we output the forces or not
964 call td_write_init(td%write_handler, namespace, space, outp, gr, st, hm, ions, ext_partners, &
965 ks, td%ions_dyn%is_active(), &
966 list_has_gauge_field(ext_partners), hm%kick, td%iter, td%max_iter, td%dt, mc)
967
968 ! Resets the nondipole integration after laser-file has been written.
969 lasers => list_get_lasers(ext_partners)
970 if(associated(lasers)) then
971 if (lasers_with_nondipole_field(lasers)) then
972 ndinitial(1:space%dim)=m_zero
973 call lasers_set_nondipole_parameters(lasers,ndinitial,m_zero)
974 end if
975 end if
976 nullify(lasers)
977
978 call td_init_ions_and_forces(td, namespace, space, gr, ions, ext_partners, st, ks, hm, outp)
979
980 if (td%scissor > m_epsilon) then
981 call scissor_init(hm%scissor, namespace, space, st, gr, hm%kpoints, hm%phase, td%scissor, mc)
982 end if
983
984 if (td%iter == 0) call td_run_zero_iter(td, namespace, space, gr, ions, st, ks, hm, ext_partners, outp, mc)
985
986 gfield => list_get_gauge_field(ext_partners)
987 if(associated(gfield)) then
988 if (gauge_field_is_propagated(gfield)) then
989 if(abs(ks%xc%kernel_lrc_alpha) > m_epsilon) then
990 call gauge_field_get_force(gfield, gr, st%d%spin_channels, st%current, &
991 ks%xc%kernel_lrc_alpha, ks%xc%kernel_proca_a_zero, ks%xc%kernel_proca_a_one)
992 call messages_experimental('TD-LRC kernel')
993 else
994 call gauge_field_get_force(gfield, gr, st%d%spin_channels, st%current)
995 endif
996 endif
997 end if
998
999 !call td_check_trotter(td, sys, h)
1000 td%iter = td%iter + 1
1001
1002 call td%restart_dump%init(namespace, restart_td, restart_type_dump, mc, ierr, mesh=gr)
1003 if (td%ions_dyn%is_active() .and. td%recalculate_gs) then
1004 ! We will also use the TD restart directory as temporary storage during the time propagation
1005 call td%restart_load%init(namespace, restart_td, restart_type_load, mc, ierr, mesh=gr)
1006 end if
1007
1008 call messages_print_with_emphasis(msg="Time-Dependent Simulation", namespace=namespace)
1009 call td_print_header(namespace)
1010
1011 if (td%pesv%calc_spm .or. td%pesv%calc_mask .and. from_scratch) then
1012 call pes_init_write(td%pesv,gr,st, namespace)
1013 end if
1014
1015 if (st%pack_states .and. hm%apply_packed()) call st%pack()
1016
1018 end subroutine td_init_with_wavefunctions
1019
1020 ! ---------------------------------------------------------
1021 subroutine td_init_ions_and_forces(td, namespace, space, gr, ions, ext_partners, st, ks, hm, outp)
1022 type(td_t), intent(inout) :: td
1023 type(namespace_t), intent(in) :: namespace
1024 type(electron_space_t), intent(in) :: space
1025 type(grid_t), intent(inout) :: gr
1026 type(ions_t), intent(inout) :: ions
1027 type(partner_list_t), intent(in) :: ext_partners
1028 type(states_elec_t), target, intent(inout) :: st
1029 type(v_ks_t), intent(inout) :: ks
1030 type(hamiltonian_elec_t), intent(inout) :: hm
1031 type(output_t), intent(inout) :: outp
1032
1033 push_sub(td_init_ions_and_forces)
1034
1035 ! Calculate initial forces and kinetic energy
1036 if (td%ions_dyn%ions_move()) then
1037 if (td%iter > 0) then
1038 call td_read_coordinates(td, namespace, ions)
1039 if (ion_dynamics_drive_ions(td%ions_dyn)) then
1040 call ion_dynamics_propagate(td%ions_dyn, ions, td%iter*td%dt, td%dt, namespace)
1041 end if
1042 call hamiltonian_elec_epot_generate(hm, namespace, space, gr, ions, ext_partners, st, time = td%iter*td%dt)
1043 ! recompute potential because the ions have moved
1044 call v_ks_calc(ks, namespace, space, hm, st, ions, ext_partners, calc_eigenval=.true., time = td%iter*td%dt)
1045 end if
1046
1047 call forces_calculate(gr, namespace, ions, hm, ext_partners, st, ks, t = td%iter*td%dt, dt = td%dt)
1048
1049 call ions%update_kinetic_energy()
1050 else
1051 if (outp%what(option__output__forces) .or. td%write_handler%out(out_separate_forces)%write) then
1052 call forces_calculate(gr, namespace, ions, hm, ext_partners, st, ks, t = td%iter*td%dt, dt = td%dt)
1053 end if
1054 end if
1055
1056 if (outp%what(option__output__stress) .or. td%ions_dyn%cell_relax()) then
1057 call stress_calculate(namespace, gr, hm, st, ions, ks, ext_partners)
1058 if (td%ions_dyn%cell_relax()) then
1059 call td%ions_dyn%update_stress(ions%space, st%stress_tensors%total, ions%latt%rlattice, ions%latt%rcell_volume)
1060 end if
1061 end if
1062
1064 end subroutine td_init_ions_and_forces
1065
1066 ! ---------------------------------------------------------
1067 subroutine td_load_restart_from_td(td, namespace, space, mc, gr, ext_partners, st, ks, hm, from_scratch)
1068 type(td_t), intent(inout) :: td
1069 type(namespace_t), intent(in) :: namespace
1070 class(space_t), intent(in) :: space
1071 type(multicomm_t), intent(in) :: mc
1072 type(grid_t), intent(inout) :: gr
1073 type(partner_list_t), intent(in) :: ext_partners
1074 type(states_elec_t), target, intent(inout) :: st
1075 type(v_ks_t), intent(inout) :: ks
1076 type(hamiltonian_elec_t), intent(inout) :: hm
1077 logical, intent(inout) :: from_scratch
1078
1079 integer :: ierr
1080 type(restart_t) :: restart
1081
1082 push_sub(td_load_restart_from_td)
1083
1084 !We redistribute the states before the restarting
1085 if (td%freeze_orbitals > 0) then
1086 call states_elec_freeze_redistribute_states(st, namespace, gr, mc, td%freeze_orbitals)
1087 end if
1088
1089 call restart%init(namespace, restart_td, restart_type_load, mc, ierr, mesh=gr)
1090 if (ierr == 0) then
1091 call td_load(restart, namespace, space, gr, st, hm, ext_partners, td, ks, ierr)
1092 end if
1093 call restart%end()
1094 if (ierr /= 0) then
1095 from_scratch = .true.
1096 td%iter = 0
1097 end if
1098
1100 end subroutine td_load_restart_from_td
1101
1102 ! ---------------------------------------------------------
1103 subroutine td_load_restart_from_gs(td, namespace, space, mc, gr, ext_partners, st, ks, hm)
1104 type(td_t), intent(inout) :: td
1105 type(namespace_t), intent(in) :: namespace
1106 class(space_t), intent(in) :: space
1107 type(multicomm_t), intent(in) :: mc
1108 type(grid_t), intent(inout) :: gr
1109 type(partner_list_t), intent(in) :: ext_partners
1110 type(states_elec_t), target, intent(inout) :: st
1111 type(v_ks_t), intent(inout) :: ks
1112 type(hamiltonian_elec_t), intent(inout) :: hm
1113
1114 integer :: ierr
1115 type(restart_t) :: restart
1117 push_sub(td_load_restart_from_gs)
1118
1119 call restart%init(namespace, restart_gs, restart_type_load, mc, ierr, mesh=gr, exact=.true.)
1120
1121 if (.not. st%only_userdef_istates) then
1122 if (ierr == 0) then
1123 call states_elec_load(restart, namespace, space, st, gr, hm%kpoints, ierr, label = ": gs")
1124 end if
1125 if (ierr /= 0) then
1126 message(1) = 'Unable to read ground-state wavefunctions.'
1127 call messages_fatal(1, namespace=namespace)
1128 end if
1129 end if
1130
1131 ! check if we should deploy user-defined wavefunctions.
1132 ! according to the settings in the input file the routine
1133 ! overwrites orbitals that were read from restart/gs
1134 if (parse_is_defined(namespace, 'UserDefinedStates')) then
1135 call states_elec_read_user_def_orbitals(gr, namespace, space, st)
1136 end if
1137
1138 call states_elec_transform(st, namespace, space, restart, gr, hm%kpoints)
1139 call restart%end()
1140
1142 end subroutine td_load_restart_from_gs
1143
1144 ! ---------------------------------------------------------
1145 subroutine td_run_zero_iter(td, namespace, space, gr, ions, st, ks, hm, ext_partners, outp, mc)
1146 type(td_t), intent(inout) :: td
1147 type(namespace_t), intent(in) :: namespace
1148 type(electron_space_t), intent(in) :: space
1149 type(grid_t), intent(inout) :: gr
1150 type(ions_t), intent(inout) :: ions
1151 type(states_elec_t), intent(inout) :: st
1152 type(v_ks_t), intent(inout) :: ks
1153 type(hamiltonian_elec_t), intent(inout) :: hm
1154 type(partner_list_t), intent(in) :: ext_partners
1155 type(output_t), intent(in) :: outp
1156 type(multicomm_t), intent(in) :: mc
1157
1158 push_sub(td_run_zero_iter)
1159
1160 call td_write_iter(td%write_handler, namespace, space, outp, gr, st, hm, ions, ext_partners, &
1161 hm%kick, ks, td%dt, 0, mc, td%recalculate_gs)
1163 ! I apply the delta electric field *after* td_write_iter, otherwise the
1164 ! dipole matrix elements in write_proj are wrong
1165 if (abs(hm%kick%time) <= m_epsilon) then
1166 if (.not. hm%pcm%localf) then
1167 call kick_apply(space, gr, st, td%ions_dyn, ions, hm%kick, hm%psolver, hm%kpoints)
1168 else
1169 call kick_apply(space, gr, st, td%ions_dyn, ions, hm%kick, hm%psolver, hm%kpoints, pcm = hm%pcm)
1170 end if
1171 call td_write_kick(outp, namespace, space, gr, hm%kick, ions, 0)
1172
1173 !We activate the sprial BC only after the kick
1174 if (gr%der%boundaries%spiralBC) then
1175 gr%der%boundaries%spiral = .true.
1176 end if
1177 end if
1178 call hm%ks_pot%run_zero_iter(td%tr%vks_old)
1179
1180 if (any(outp%output_interval > 0)) then
1181 call td_write_data(td%write_handler)
1182 call td_write_output(namespace, space, gr, st, hm, ks, outp, ions, ext_partners, 0)
1183 end if
1184
1185 pop_sub(td_run_zero_iter)
1186 end subroutine td_run_zero_iter
1187
1188
1189 ! ---------------------------------------------------------
1191 subroutine td_read_coordinates(td, namespace, ions)
1192 type(td_t), intent(in) :: td
1193 type(namespace_t), intent(in) :: namespace
1194 type(ions_t), intent(inout) :: ions
1195
1196 integer :: iatom, iter, iunit
1197
1199
1200 iunit = io_open('td.general/coordinates', namespace, action='read', status='old', die=.false.)
1201 if (iunit == -1) then
1202 message(1) = "Could not open file '"//trim(io_workpath('td.general/coordinates', namespace))//"'."
1203 message(2) = "Starting simulation from initial geometry."
1204 call messages_warning(2, namespace=namespace)
1205 pop_sub(td_read_coordinates)
1206 return
1207 end if
1208
1209 call io_skip_header(iunit)
1210 do iter = 0, td%iter - 1
1211 read(iunit, *) ! skip previous iterations... sorry, but no portable seek in Fortran
1212 end do
1213 read(iunit, '(32x)', advance='no') ! skip the time index.
1214
1215 do iatom = 1, ions%natoms
1216 read(iunit, '(3es24.16)', advance='no') ions%pos(:, iatom)
1217 ions%pos(:, iatom) = units_to_atomic(units_out%length, ions%pos(:, iatom))
1218 end do
1219 do iatom = 1, ions%natoms
1220 read(iunit, '(3es24.16)', advance='no') ions%vel(:, iatom)
1221 ions%vel(:, iatom) = units_to_atomic(units_out%velocity, ions%vel(:, iatom))
1222 end do
1223 do iatom = 1, ions%natoms
1224 read(iunit, '(3es24.16)', advance='no') ions%tot_force(:, iatom)
1225 ions%tot_force(:, iatom) = units_to_atomic(units_out%force, ions%tot_force(:, iatom))
1226 end do
1227
1228 call io_close(iunit)
1229
1230 pop_sub(td_read_coordinates)
1231 end subroutine td_read_coordinates
1232
1233 ! ---------------------------------------------------------
1234 subroutine td_dump(td, namespace, space, gr, st, hm, ks, ext_partners, iter, ierr)
1235 type(td_t), intent(in) :: td
1236 type(namespace_t), intent(in) :: namespace
1237 class(space_t), intent(in) :: space
1238 type(grid_t), intent(in) :: gr
1239 type(states_elec_t), intent(in) :: st
1240 type(hamiltonian_elec_t), intent(in) :: hm
1241 type(v_ks_t), intent(in) :: ks
1242 type(partner_list_t), intent(in) :: ext_partners
1243 integer, intent(in) :: iter
1244 integer, intent(out) :: ierr
1245
1246 type(gauge_field_t), pointer :: gfield
1247 integer :: err, err2
1248
1249 push_sub(td_dump)
1250
1251 ierr = 0
1252
1253 if (td%restart_dump%skip()) then
1254 pop_sub(td_dump)
1255 return
1256 end if
1257
1258 message(1) = "Debug: Writing td restart."
1259 call messages_info(1, namespace=namespace, debug_only=.true.)
1260
1261 ! first write resume file
1262 call states_elec_dump(td%restart_dump, space, st, gr, hm%kpoints, err, iter=iter)
1263 if (err /= 0) ierr = ierr + 1
1264
1265 call states_elec_dump_rho(td%restart_dump, space, st, gr, ierr, iter=iter)
1266 if (err /= 0) ierr = ierr + 1
1267
1268 if (hm%lda_u_level /= dft_u_none) then
1269 call lda_u_dump(td%restart_dump, namespace, hm%lda_u, st, gr, ierr)
1270 if (err /= 0) ierr = ierr + 1
1271 end if
1272
1273 call potential_interpolation_dump(td%tr%vks_old, space, td%restart_dump, gr, st%d%nspin, err2)
1274 if (err2 /= 0) ierr = ierr + 2
1275
1276 call pes_dump(td%pesv, namespace, td%restart_dump, st, gr, err)
1277 if (err /= 0) ierr = ierr + 4
1278
1279 ! Gauge field restart
1280 gfield => list_get_gauge_field(ext_partners)
1281 if(associated(gfield)) then
1282 call gauge_field_dump(td%restart_dump, gfield, ierr)
1283 end if
1284
1285 if (gr%der%boundaries%spiralBC) then
1286 call states_elec_dump_spin(td%restart_dump, st, err)
1287 if (err /= 0) ierr = ierr + 8
1288 end if
1289
1290 if (ks%has_photons) then
1291 call mf_photons_dump(td%restart_dump, ks%pt_mx, gr, td%dt, ks%pt, err)
1292 if (err /= 0) ierr = ierr + 16
1293 end if
1294
1295 if (ks%xc_photon /= 0) then
1296 ! photon-free mean field
1297 call ks%xc_photons%mf_dump(td%restart_dump, err)
1298 if (err /= 0) ierr = ierr + 32
1299 end if
1300
1301 if (allocated(st%frozen_rho)) then
1302 call states_elec_dump_frozen(td%restart_dump, space, st, gr, ierr)
1303 end if
1304 if (err /= 0) ierr = ierr + 64
1305
1306 if (td%ions_dyn%ions_move() .or. td%ions_dyn%cell_relax()) then
1307 call ion_dynamics_dump(td%ions_dyn, td%restart_dump, err)
1308 end if
1309 if (err /= 0) ierr = ierr + 128
1310
1311 message(1) = "Debug: Writing td restart done."
1312 call messages_info(1, namespace=namespace, debug_only=.true.)
1313
1314 pop_sub(td_dump)
1315 end subroutine td_dump
1316
1317 ! ---------------------------------------------------------
1318 subroutine td_load(restart, namespace, space, gr, st, hm, ext_partners, td, ks, ierr)
1319 type(restart_t), intent(in) :: restart
1320 type(namespace_t), intent(in) :: namespace
1321 class(space_t), intent(in) :: space
1322 type(grid_t), intent(in) :: gr
1323 type(states_elec_t), intent(inout) :: st
1324 type(hamiltonian_elec_t), intent(inout) :: hm
1325 type(partner_list_t), intent(in) :: ext_partners
1326 type(td_t), intent(inout) :: td
1327 type(v_ks_t), intent(inout) :: ks
1328 integer, intent(out) :: ierr
1330 integer :: err, err2
1331 type(gauge_field_t), pointer :: gfield
1332 push_sub(td_load)
1333
1334 ierr = 0
1335
1336 if (restart%skip()) then
1337 ierr = -1
1338 pop_sub(td_load)
1339 return
1340 end if
1341
1342 message(1) = "Debug: Reading td restart."
1343 call messages_info(1, namespace=namespace, debug_only=.true.)
1344
1345 ! Read states
1346 call states_elec_load(restart, namespace, space, st, gr, hm%kpoints, err, iter=td%iter, label = ": td")
1347 if (err /= 0) then
1348 ierr = ierr + 1
1349 end if
1350
1351 ! read potential from previous interactions
1352 call potential_interpolation_load(td%tr%vks_old, namespace, space, restart, gr, st%d%nspin, err2)
1353 if (err2 /= 0) ierr = ierr + 2
1354
1355 if (hm%lda_u_level /= dft_u_none) then
1356 call lda_u_load(restart, hm%lda_u, st, hm%energy%dft_u, ierr)
1357 if (err /= 0) ierr = ierr + 1
1358 end if
1359
1360
1361 ! read PES restart
1362 if (td%pesv%calc_spm .or. td%pesv%calc_mask .or. td%pesv%calc_flux) then
1363 call pes_load(td%pesv, namespace, restart, st, err)
1364 if (err /= 0) ierr = ierr + 4
1365 end if
1366
1367 ! Gauge field restart
1368 gfield => list_get_gauge_field(ext_partners)
1369 if (associated(gfield)) then
1370 call gauge_field_load(restart, gfield, err)
1371 if (err /= 0) then
1372 ierr = ierr + 8
1373 else
1374 call hm%update(gr, namespace, space, ext_partners, time = td%dt*td%iter)
1375 end if
1376 end if
1377
1378 ! add photon restart
1379 if (ks%has_photons) then
1380 call mf_photons_load(restart, ks%pt_mx, gr, err)
1381 end if
1382 if (err /= 0) ierr = ierr + 16
1383
1384 if (ks%xc_photon /= 0) then
1385 call ks%xc_photons%mf_load(restart, space, err)
1386 end if
1387 if (err /= 0) ierr = ierr + 32
1388
1389 if (gr%der%boundaries%spiralBC) then
1390 call states_elec_load_spin(restart, st, err)
1391 !To ensure back compatibility, if the file is not present, we use the
1392 !current states to get the spins
1393 if (err /= 0) call states_elec_fermi(st, namespace, gr)
1394 end if
1395
1396 if (td%ions_dyn%is_active()) then
1397 call ion_dynamics_load(td%ions_dyn, restart, err)
1398 end if
1399 if (err /= 0) ierr = ierr + 64
1400
1401 message(1) = "Debug: Reading td restart done."
1402 call messages_info(1, namespace=namespace, debug_only=.true.)
1403
1404 pop_sub(td_load)
1405 end subroutine td_load
1406 ! ---------------------------------------------------------
1407 subroutine td_load_frozen(namespace, restart, space, mesh, st, hm, ierr)
1408 type(namespace_t), intent(in) :: namespace
1409 type(restart_t), intent(in) :: restart
1410 class(space_t), intent(in) :: space
1411 class(mesh_t), intent(in) :: mesh
1412 type(states_elec_t), intent(inout) :: st
1413 type(hamiltonian_elec_t), intent(inout) :: hm
1414 integer, intent(out) :: ierr
1415
1416 push_sub(td_load_frozen)
1417
1418 ierr = 0
1419
1420 if (restart%skip()) then
1421 ierr = -1
1422 pop_sub(td_load_frozen)
1423 return
1424 end if
1425
1426 message(1) = "Debug: Reading td frozen restart."
1427 call messages_info(1, namespace=namespace, debug_only=.true.)
1428
1429 safe_allocate(st%frozen_rho(1:mesh%np, 1:st%d%nspin))
1430 if (family_is_mgga(hm%xc%family)) then
1431 safe_allocate(st%frozen_tau(1:mesh%np, 1:st%d%nspin))
1432 safe_allocate(st%frozen_gdens(1:mesh%np, 1:space%dim, 1:st%d%nspin))
1433 safe_allocate(st%frozen_ldens(1:mesh%np, 1:st%d%nspin))
1434 end if
1435
1436 call states_elec_load_frozen(restart, space, st, mesh, ierr)
1437
1438 message(1) = "Debug: Reading td frozen restart done."
1439 call messages_info(1, namespace=namespace, debug_only=.true.)
1440
1441 pop_sub(td_load_frozen)
1442 end subroutine td_load_frozen
1443
1444 ! ---------------------------------------------------------
1445 logical function td_get_from_scratch(td)
1446 type(td_t), intent(in) :: td
1447
1448 push_sub(td_get_from_scratch)
1449
1450 td_get_from_scratch = td%from_scratch
1451
1452 pop_sub(td_get_from_scratch)
1453 end function td_get_from_scratch
1454
1455 ! ---------------------------------------------------------
1456 subroutine td_set_from_scratch(td, from_scratch)
1457 type(td_t), intent(inout) :: td
1458 logical, intent(in) :: from_scratch
1459
1460 push_sub(td_set_from_scratch)
1461
1462 td%from_scratch = from_scratch
1463
1464 pop_sub(td_set_from_scratch)
1465 end subroutine td_set_from_scratch
1466end module td_oct_m
1467
1468!! Local Variables:
1469!! mode: f90
1470!! coding: utf-8
1471!! 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_alpha, proca_a_zero, proca_a_one)
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, space, restart, mesh, nspin, err2)
subroutine, public potential_interpolation_dump(potential_interpolation, space, 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:255
subroutine, public scf_end(scf)
Definition: scf.F90:525
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, space, 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:787
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:1199
subroutine td_print_message(td, namespace, ions, hm, iter, scsteps, etime)
Definition: td.F90:851
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:1541
subroutine td_load_frozen(namespace, restart, space, mesh, st, hm, ierr)
Definition: td.F90:1503
subroutine td_print_header(namespace)
Definition: td.F90:772
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:1414
subroutine, public td_set_from_scratch(td, from_scratch)
Definition: td.F90:1552
subroutine, public td_dump(td, namespace, space, gr, st, hm, ks, ext_partners, iter, ierr)
Definition: td.F90:1330
subroutine td_read_coordinates(td, namespace, ions)
reads the pos and vel from coordinates file
Definition: td.F90:1287
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:1241
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:872
subroutine, public td_init_with_wavefunctions(td, namespace, space, mc, gr, ions, ext_partners, st, ks, hm, outp, from_scratch)
Definition: td.F90:883
subroutine, public td_load_restart_from_td(td, namespace, space, mc, gr, ext_partners, st, ks, hm, from_scratch)
Definition: td.F90:1163
subroutine td_init_ions_and_forces(td, namespace, space, gr, ions, ext_partners, st, ks, hm, outp)
Definition: td.F90:1117
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:628
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:641
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)