29 use,
intrinsic :: iso_fortran_env
74 type(space_t),
public :: space
75 real(real64),
allocatable :: vecpot(:)
76 real(real64),
allocatable :: vecpot_vel(:)
77 real(real64),
allocatable :: vecpot_acc(:)
78 real(real64),
allocatable :: vecpot_kick(:)
79 real(real64),
allocatable :: force(:)
81 logical :: with_gauge_field = .false.
83 real(real64) :: kicktime
85 real(real64) :: volume
102 class(gauge_field_t),
pointer :: this
103 type(namespace_t),
intent(in) :: namespace
104 real(real64),
intent(in) :: volume
108 logical :: default_gauge_field_propagate
114 this%namespace =
namespace_t(
"GaugeField", parent=namespace)
118 this%space =
space_t(namespace)
120 this%with_gauge_field = .false.
124 safe_allocate(this%vecpot(1:this%space%dim))
125 safe_allocate(this%vecpot_vel(1:this%space%dim))
126 safe_allocate(this%vecpot_acc(1:this%space%dim))
127 safe_allocate(this%vecpot_kick(1:this%space%dim))
128 safe_allocate(this%force(1:this%space%dim))
149 call parse_variable(namespace,
'GaugeFieldDynamics', option__gaugefielddynamics__polarization, this%dynamics)
172 default_gauge_field_propagate = .false.
175 do ii = 1, this%space%dim
180 if (norm2(this%vecpot_kick(1:this%space%dim)) >
m_epsilon)
then
181 default_gauge_field_propagate = .
true.
185 if (.not. this%space%is_periodic())
then
186 message(1) =
"GaugeVectorField is intended for periodic systems."
199 call parse_variable(namespace,
'GaugeFieldPropagate', default_gauge_field_propagate, this%with_gauge_field)
212 if (abs(this%kicktime) <=
m_epsilon)
then
213 this%vecpot(1:this%space%dim) = this%vecpot_kick(1:this%space%dim)
228 if (kpoints%use_symmetries)
then
232 message(1) =
"The GaugeVectorField breaks (at least) one of the symmetries used to reduce the k-points."
233 message(2) =
"Set SymmetryBreakDir equal to GaugeVectorField."
249 this%with_gauge_field = .false.
250 safe_deallocate_a(this%vecpot)
251 safe_deallocate_a(this%vecpot_vel)
252 safe_deallocate_a(this%vecpot_acc)
253 safe_deallocate_a(this%vecpot_kick)
254 safe_deallocate_a(this%force)
264 is_propagated = this%with_gauge_field
271 is_used = this%with_gauge_field .or. (norm2(this%vecpot_kick(1:this%space%dim)) > m_epsilon)
279 real(real64),
intent(in) :: vec_pot(:)
282 this%vecpot(1:this%space%dim) = vec_pot(1:this%space%dim)
291 real(real64),
intent(in) :: vec_pot_vel(:)
294 this%vecpot_vel(1:this%space%dim) = vec_pot_vel(1:this%space%dim)
303 real(real64),
intent(out) :: vec_pot(:)
306 vec_pot(1:this%space%dim) = this%vecpot(1:this%space%dim)
315 real(real64),
intent(out) :: vec_pot_vel(:)
318 vec_pot_vel(1:this%space%dim) = this%vecpot_vel(1:this%space%dim)
327 real(real64),
intent(out) :: vec_pot_acc(:)
330 vec_pot_acc(1:this%space%dim) = this%vecpot_acc(1:this%space%dim)
338 real(real64),
intent(in) :: dt
339 real(real64),
intent(in) :: time
341 logical,
save :: warning_shown = .false.
346 this%vecpot_acc(1:this%space%dim) = this%force(1:this%space%dim)
349 if (this%kicktime > m_zero .and. time-dt <= this%kicktime .and. time > this%kicktime)
then
350 this%vecpot(1:this%space%dim) = this%vecpot(1:this%space%dim) + this%vecpot_kick(1:this%space%dim)
351 call messages_write(
' ---------------- Applying gauge kick ----------------')
352 call messages_info(namespace=this%namespace)
355 this%vecpot(1:this%space%dim) = this%vecpot(1:this%space%dim) + dt * this%vecpot_vel(1:this%space%dim) + &
356 m_half * dt**2 * this%force(1:this%space%dim)
359 do idim = 1, this%space%dim
360 if (.not. warning_shown .and. abs(this%vecpot_kick(idim)) > m_epsilon .and. &
361 abs(this%vecpot(idim))> abs(this%vecpot_kick(idim))*1.01 .and. .not. this%kicktime > m_zero)
then
365 write(message(1),
'(a)')
'It seems that the gauge-field might be diverging. You should probably check'
366 write(message(2),
'(a)')
'the simulation parameters, in particular the number of k-points.'
367 call messages_warning(2, namespace=this%namespace)
376 real(real64),
intent(in) :: qtot
380 this%wp2 = m_four*m_pi*qtot/this%volume
382 write (message(1),
'(a,f12.6,a)')
"Info: Electron-gas plasmon frequency", &
383 units_from_atomic(units_out%energy,
sqrt(this%wp2)),
" ["//trim(units_abbrev(units_out%energy))//
"]"
384 call messages_info(1, namespace=this%namespace)
395 if(
allocated(this%vecpot_vel))
then
396 energy = this%volume / (8.0_real64 * m_pi * p_c**2) * sum(this%vecpot_vel(1:this%space%dim)**2)
407 type(restart_t),
intent(in) :: restart
409 integer,
intent(out) :: ierr
412 real(real64),
allocatable :: vecpot(:,:)
418 if (restart%skip())
then
423 message(1) =
"Debug: Writing gauge field restart."
424 call messages_info(1, namespace=gfield%namespace, debug_only=.
true.)
426 safe_allocate(vecpot(1:gfield%space%dim, 1:2))
431 call restart%write_binary(
"gauge_field", 2*gfield%space%dim, vecpot, err)
432 safe_deallocate_a(vecpot)
433 if (err /= 0) ierr = ierr + 1
435 message(1) =
"Debug: Writing gauge field restart done."
436 call messages_info(1, namespace=gfield%namespace, debug_only=.
true.)
444 type(restart_t),
intent(in) :: restart
446 integer,
intent(out) :: ierr
449 real(real64),
allocatable :: vecpot(:,:)
455 if (restart%skip())
then
461 message(1) =
"Debug: Reading gauge field restart."
462 call messages_info(1, namespace=gfield%namespace, debug_only=.
true.)
464 safe_allocate(vecpot(1:gfield%space%dim, 1:2))
465 call restart%read_binary(
"gauge_field", 2*gfield%space%dim, vecpot, err)
466 if (err /= 0) ierr = ierr + 1
470 safe_deallocate_a(vecpot)
472 message(1) =
"Debug: Reading gauge field restart done."
473 call messages_info(1, namespace=gfield%namespace, debug_only=.
true.)
482 type(grid_t),
intent(in) :: gr
483 integer,
intent(in) :: spin_channels
484 real(real64),
intent(in) :: current(:,:,:)
485 type(xc_lrc_t),
intent(in) :: lrc
487 integer :: idir, ispin
491 select case (this%dynamics)
492 case (option__gaugefielddynamics__none)
493 this%force(1:this%space%dim) = m_zero
495 case (option__gaugefielddynamics__polarization)
496 do idir = 1, this%space%periodic_dim
497 this%force(idir) = m_zero
498 do ispin = 1, spin_channels
503 if(abs(lrc%alpha) > m_epsilon)
then
504 this%force(idir) = this%force(idir) + &
505 lrc%alpha*p_c/this%volume*dmf_integrate(gr, current(:, idir, ispin)) + &
506 lrc%alpha*lrc%proca_a_zero*this%vecpot_kick(idir) - &
507 lrc%alpha*lrc%proca_a_zero*this%vecpot(idir) - lrc%alpha*lrc%proca_a_one*this%vecpot_vel(idir)
509 this%force(idir) = this%force(idir) - &
510 m_four*m_pi*p_c/this%volume*dmf_integrate(gr, current(:, idir, ispin))
526 type(algorithmic_operation_t),
intent(in) :: operation
527 real(real64),
intent(in) :: dt
528 real(real64),
intent(in) :: time
532 select case (operation%id)
537 case (verlet_update_pos)
539 case (verlet_compute_acc)
542 case (verlet_compute_vel)
543 this%vecpot_vel(1:this%space%dim) = this%vecpot_vel(1:this%space%dim) + &
544 m_half * dt * (this%vecpot_acc(1:this%space%dim) + this%force(1:this%space%dim))
546 message(1) =
"Unsupported TD operation."
547 call messages_fatal(1, namespace=this%namespace)
557 type(c_ptr),
intent(inout) :: out_gauge
558 integer,
intent(in) :: iter
561 character(len=50) :: aux
562 real(real64) :: temp(this%space%dim)
564 if (.not. mpi_world%is_root())
return
569 call write_iter_clear(out_gauge)
570 call write_iter_string(out_gauge,
'################################################################################')
571 call write_iter_nl(out_gauge)
572 call write_iter_string(out_gauge,
'# HEADER')
573 call write_iter_nl(out_gauge)
576 call write_iter_header_start(out_gauge)
578 do idir = 1, this%space%dim
579 write(aux,
'(a2,i1,a1)')
'A(', idir,
')'
580 call write_iter_header(out_gauge, aux)
582 do idir = 1, this%space%dim
583 write(aux,
'(a6,i1,a1)')
'dA/dt(', idir,
')'
584 call write_iter_header(out_gauge, aux)
586 do idir = 1, this%space%dim
587 write(aux,
'(a10,i1,a1)')
'd^2A/dt^2(', idir,
')'
588 call write_iter_header(out_gauge, aux)
590 call write_iter_nl(out_gauge)
601 call write_iter_string(out_gauge,
'################################################################################')
602 call write_iter_nl(out_gauge)
606 call write_iter_start(out_gauge)
608 do idir = 1, this%space%dim
609 temp(idir) = units_from_atomic(units_out%energy, this%vecpot(idir))
611 call write_iter_double(out_gauge, temp, this%space%dim)
613 do idir = 1, this%space%dim
614 temp(idir) = units_from_atomic(units_out%energy / units_out%time, this%vecpot_vel(idir))
616 call write_iter_double(out_gauge, temp, this%space%dim)
618 do idir = 1, this%space%dim
619 temp(idir) = units_from_atomic(units_out%energy / units_out%time**2, this%vecpot_acc(idir))
621 call write_iter_double(out_gauge, temp, this%space%dim)
623 call write_iter_nl(out_gauge)
642 class(interaction_surrogate_t),
intent(inout) :: interaction
646 select type (interaction)
648 message(1) =
"Unsupported interaction."
649 call messages_fatal(1)
658 class(interaction_surrogate_t),
intent(inout) :: interaction
662 select type (interaction)
664 message(1) =
"Unsupported interaction."
665 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_get_force(this, gr, spin_channels, current, lrc)
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)
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