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