29 use,
intrinsic :: iso_fortran_env
73 type(space_t),
public :: space
74 real(real64),
allocatable :: vecpot(:)
75 real(real64),
allocatable :: vecpot_vel(:)
76 real(real64),
allocatable :: vecpot_acc(:)
77 real(real64),
allocatable :: vecpot_kick(:)
78 real(real64),
allocatable :: force(:)
80 logical :: with_gauge_field = .false.
82 real(real64) :: kicktime
84 real(real64) :: volume
101 class(gauge_field_t),
pointer :: this
102 type(namespace_t),
intent(in) :: namespace
103 real(real64),
intent(in) :: volume
107 logical :: default_gauge_field_propagate
113 this%namespace =
namespace_t(
"GaugeField", parent=namespace)
119 this%with_gauge_field = .false.
123 safe_allocate(this%vecpot(1:this%space%dim))
124 safe_allocate(this%vecpot_vel(1:this%space%dim))
125 safe_allocate(this%vecpot_acc(1:this%space%dim))
126 safe_allocate(this%vecpot_kick(1:this%space%dim))
127 safe_allocate(this%force(1:this%space%dim))
148 call parse_variable(namespace,
'GaugeFieldDynamics', option__gaugefielddynamics__polarization, this%dynamics)
171 default_gauge_field_propagate = .false.
174 do ii = 1, this%space%dim
179 if (norm2(this%vecpot_kick(1:this%space%dim)) >
m_epsilon)
then
180 default_gauge_field_propagate = .
true.
184 if (.not. this%space%is_periodic())
then
185 message(1) =
"GaugeVectorField is intended for periodic systems."
198 call parse_variable(namespace,
'GaugeFieldPropagate', default_gauge_field_propagate, this%with_gauge_field)
211 if (abs(this%kicktime) <=
m_epsilon)
then
212 this%vecpot(1:this%space%dim) = this%vecpot_kick(1:this%space%dim)
227 if (kpoints%use_symmetries)
then
231 message(1) =
"The GaugeVectorField breaks (at least) one of the symmetries used to reduce the k-points."
232 message(2) =
"Set SymmetryBreakDir equal to GaugeVectorField."
248 this%with_gauge_field = .false.
249 safe_deallocate_a(this%vecpot)
250 safe_deallocate_a(this%vecpot_vel)
251 safe_deallocate_a(this%vecpot_acc)
252 safe_deallocate_a(this%vecpot_kick)
253 safe_deallocate_a(this%force)
263 is_propagated = this%with_gauge_field
270 is_used = this%with_gauge_field .or. (norm2(this%vecpot_kick(1:this%space%dim)) > m_epsilon)
278 real(real64),
intent(in) :: vec_pot(:)
281 this%vecpot(1:this%space%dim) = vec_pot(1:this%space%dim)
290 real(real64),
intent(in) :: vec_pot_vel(:)
293 this%vecpot_vel(1:this%space%dim) = vec_pot_vel(1:this%space%dim)
302 real(real64),
intent(out) :: vec_pot(:)
305 vec_pot(1:this%space%dim) = this%vecpot(1:this%space%dim)
314 real(real64),
intent(out) :: vec_pot_vel(:)
317 vec_pot_vel(1:this%space%dim) = this%vecpot_vel(1:this%space%dim)
326 real(real64),
intent(out) :: vec_pot_acc(:)
329 vec_pot_acc(1:this%space%dim) = this%vecpot_acc(1:this%space%dim)
337 real(real64),
intent(in) :: dt
338 real(real64),
intent(in) :: time
340 logical,
save :: warning_shown = .false.
345 this%vecpot_acc(1:this%space%dim) = this%force(1:this%space%dim)
348 if (this%kicktime > m_zero .and. time-dt <= this%kicktime .and. time > this%kicktime)
then
349 this%vecpot(1:this%space%dim) = this%vecpot(1:this%space%dim) + this%vecpot_kick(1:this%space%dim)
350 call messages_write(
' ---------------- Applying gauge kick ----------------')
351 call messages_info(namespace=this%namespace)
354 this%vecpot(1:this%space%dim) = this%vecpot(1:this%space%dim) + dt * this%vecpot_vel(1:this%space%dim) + &
355 m_half * dt**2 * this%force(1:this%space%dim)
358 do idim = 1, this%space%dim
359 if (.not. warning_shown .and. abs(this%vecpot_kick(idim)) > m_epsilon .and. &
360 abs(this%vecpot(idim))> abs(this%vecpot_kick(idim))*1.01 .and. .not. this%kicktime > m_zero)
then
364 write(message(1),
'(a)')
'It seems that the gauge-field might be diverging. You should probably check'
365 write(message(2),
'(a)')
'the simulation parameters, in particular the number of k-points.'
366 call messages_warning(2, namespace=this%namespace)
375 real(real64),
intent(in) :: qtot
379 this%wp2 = m_four*m_pi*qtot/this%volume
381 write (message(1),
'(a,f12.6,a)')
"Info: Electron-gas plasmon frequency", &
382 units_from_atomic(units_out%energy,
sqrt(this%wp2)),
" ["//trim(units_abbrev(units_out%energy))//
"]"
383 call messages_info(1, namespace=this%namespace)
394 if(
allocated(this%vecpot_vel))
then
395 energy = this%volume / (8.0_real64 * m_pi * p_c**2) * sum(this%vecpot_vel(1:this%space%dim)**2)
406 type(restart_t),
intent(in) :: restart
408 integer,
intent(out) :: ierr
411 real(real64),
allocatable :: vecpot(:,:)
417 if (restart%skip())
then
422 message(1) =
"Debug: Writing gauge field restart."
423 call messages_info(1, namespace=gfield%namespace, debug_only=.
true.)
425 safe_allocate(vecpot(1:gfield%space%dim, 1:2))
430 call restart%write_binary(
"gauge_field", 2*gfield%space%dim, vecpot, err)
431 safe_deallocate_a(vecpot)
432 if (err /= 0) ierr = ierr + 1
434 message(1) =
"Debug: Writing gauge field restart done."
435 call messages_info(1, namespace=gfield%namespace, debug_only=.
true.)
443 type(restart_t),
intent(in) :: restart
445 integer,
intent(out) :: ierr
448 real(real64),
allocatable :: vecpot(:,:)
454 if (restart%skip())
then
460 message(1) =
"Debug: Reading gauge field restart."
461 call messages_info(1, namespace=gfield%namespace, debug_only=.
true.)
463 safe_allocate(vecpot(1:gfield%space%dim, 1:2))
464 call restart%read_binary(
"gauge_field", 2*gfield%space%dim, vecpot, err)
465 if (err /= 0) ierr = ierr + 1
469 safe_deallocate_a(vecpot)
471 message(1) =
"Debug: Reading gauge field restart done."
472 call messages_info(1, namespace=gfield%namespace, debug_only=.
true.)
479 subroutine gauge_field_get_force(this, gr, spin_channels, current, lrc_alpha, proca_a_zero, proca_a_one)
481 type(grid_t),
intent(in) :: gr
482 integer,
intent(in) :: spin_channels
483 real(real64),
intent(in) :: current(:,:,:)
484 real(real64),
optional,
intent(in) :: lrc_alpha
485 real(real64),
optional,
intent(in) :: proca_a_zero
486 real(real64),
optional,
intent(in) :: proca_a_one
488 integer :: idir, ispin
489 real(real64) :: lrc_alpha_
490 real(real64) :: proca_a_zero_
491 real(real64) :: proca_a_one_
495 lrc_alpha_ = optional_default(lrc_alpha, m_zero)
496 proca_a_zero_ = optional_default(proca_a_zero, m_zero)
497 proca_a_one_ = optional_default(proca_a_one, m_zero)
498 select case (this%dynamics)
499 case (option__gaugefielddynamics__none)
500 this%force(1:this%space%dim) = m_zero
502 case (option__gaugefielddynamics__polarization)
503 do idir = 1, this%space%periodic_dim
504 this%force(idir) = m_zero
505 do ispin = 1, spin_channels
510 if(abs(lrc_alpha_) > m_epsilon)
then
511 this%force(idir) = this%force(idir) + &
512 lrc_alpha*p_c/this%volume*dmf_integrate(gr, current(:, idir, ispin)) + &
513 lrc_alpha*proca_a_zero_*this%vecpot_kick(idir) - &
514 lrc_alpha*proca_a_zero_*this%vecpot(idir) - lrc_alpha*proca_a_one_*this%vecpot_vel(idir)
516 this%force(idir) = this%force(idir) - &
517 m_four*m_pi*p_c/this%volume*dmf_integrate(gr, current(:, idir, ispin))
533 type(algorithmic_operation_t),
intent(in) :: operation
534 real(real64),
intent(in) :: dt
535 real(real64),
intent(in) :: time
539 select case (operation%id)
544 case (verlet_update_pos)
546 case (verlet_compute_acc)
549 case (verlet_compute_vel)
550 this%vecpot_vel(1:this%space%dim) = this%vecpot_vel(1:this%space%dim) + &
551 m_half * dt * (this%vecpot_acc(1:this%space%dim) + this%force(1:this%space%dim))
553 message(1) =
"Unsupported TD operation."
554 call messages_fatal(1, namespace=this%namespace)
564 type(c_ptr),
intent(inout) :: out_gauge
565 integer,
intent(in) :: iter
568 character(len=50) :: aux
569 real(real64) :: temp(this%space%dim)
571 if (.not. mpi_world%is_root())
return
576 call write_iter_clear(out_gauge)
577 call write_iter_string(out_gauge,
'################################################################################')
578 call write_iter_nl(out_gauge)
579 call write_iter_string(out_gauge,
'# HEADER')
580 call write_iter_nl(out_gauge)
583 call write_iter_header_start(out_gauge)
585 do idir = 1, this%space%dim
586 write(aux,
'(a2,i1,a1)')
'A(', idir,
')'
587 call write_iter_header(out_gauge, aux)
589 do idir = 1, this%space%dim
590 write(aux,
'(a6,i1,a1)')
'dA/dt(', idir,
')'
591 call write_iter_header(out_gauge, aux)
593 do idir = 1, this%space%dim
594 write(aux,
'(a10,i1,a1)')
'd^2A/dt^2(', idir,
')'
595 call write_iter_header(out_gauge, aux)
597 call write_iter_nl(out_gauge)
608 call write_iter_string(out_gauge,
'################################################################################')
609 call write_iter_nl(out_gauge)
613 call write_iter_start(out_gauge)
615 do idir = 1, this%space%dim
616 temp(idir) = units_from_atomic(units_out%energy, this%vecpot(idir))
618 call write_iter_double(out_gauge, temp, this%space%dim)
620 do idir = 1, this%space%dim
621 temp(idir) = units_from_atomic(units_out%energy / units_out%time, this%vecpot_vel(idir))
623 call write_iter_double(out_gauge, temp, this%space%dim)
625 do idir = 1, this%space%dim
626 temp(idir) = units_from_atomic(units_out%energy / units_out%time**2, this%vecpot_acc(idir))
628 call write_iter_double(out_gauge, temp, this%space%dim)
630 call write_iter_nl(out_gauge)
649 class(interaction_surrogate_t),
intent(inout) :: interaction
653 select type (interaction)
655 message(1) =
"Unsupported interaction."
656 call messages_fatal(1)
665 class(interaction_surrogate_t),
intent(inout) :: interaction
669 select type (interaction)
671 message(1) =
"Unsupported interaction."
672 call messages_fatal(1)
double sqrt(double __x) __attribute__((__nothrow__
This module implements the basic elements defining algorithms.
subroutine, public gauge_field_set_vec_pot_vel(this, vec_pot_vel)
subroutine, public gauge_field_output_write(this, out_gauge, iter)
subroutine, public gauge_field_set_vec_pot(this, vec_pot)
subroutine gauge_field_finalize(this)
subroutine, public gauge_field_get_vec_pot_acc(this, vec_pot_acc)
subroutine, public gauge_field_load(restart, gfield, ierr)
subroutine gauge_field_copy_quantities_to_interaction(partner, interaction)
subroutine, public gauge_field_end(this)
subroutine, public gauge_field_get_vec_pot(this, vec_pot)
subroutine, public gauge_field_do_algorithmic_operation(this, operation, dt, time)
subroutine, public gauge_field_check_symmetries(this, kpoints)
subroutine, public gauge_field_get_vec_pot_vel(this, vec_pot_vel)
subroutine, public gauge_field_get_force(this, gr, spin_channels, current, lrc_alpha, proca_a_zero, proca_a_one)
class(gauge_field_t) function, pointer, public gauge_field_init(namespace, volume)
subroutine, public gauge_field_dump(restart, gfield, ierr)
subroutine gauge_field_init_interaction_as_partner(partner, interaction)
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) function, public gauge_field_get_energy(this)
subroutine gauge_field_propagate(this, dt, time)
real(real64), parameter, public m_zero
real(real64), parameter, public m_epsilon
This module implements the underlying real-space grid.
This module defines classes and functions for interaction partners.
This module defines various routines, operating on mesh functions.
This module defines the meshes, which are used in Octopus.
subroutine, public messages_warning(no_lines, all_nodes, namespace)
character(len=256), dimension(max_lines), public message
to be output by fatal, warning
subroutine, public messages_fatal(no_lines, only_root_writes, namespace)
integer function, public parse_block(namespace, name, blk, check_varinfo_)
integer pure function, public symmetries_identity_index(this)
integer pure function, public symmetries_number(this)
brief This module defines the class unit_t which is used by the unit_systems_oct_m module.
This module defines the unit system, used for input and output.
Explicit interfaces to C functions, defined in write_iter_low.cc.
abstract class for general interaction partners