53 type(space_t) :: space
55 integer,
public ::
type
57 character(len=1024) :: potential_formula
58 character(len=200) :: density_formula
59 character(len=MAX_PATH_LEN) :: filename
62 real(real64),
allocatable,
public :: pot(:)
64 real(real64),
allocatable,
public :: b_field(:)
66 real(real64),
allocatable,
public :: a_static(:,:)
67 real(real64),
allocatable,
public :: e_field(:)
70 real(real64),
allocatable,
public :: v_ext(:)
81 integer,
public,
parameter :: &
82 EXTERNAL_POT_USDEF = 201, & !< user-defined function for local potential
96 class(external_potential_t),
pointer :: this
97 type(namespace_t),
intent(in) :: namespace
103 this%namespace =
namespace_t(
"ExternalPotential", parent=namespace)
104 this%space =
space_t(namespace)
108 allocate(this%supported_interactions_as_partner(0))
110 call this%quantities%add(
quantity_t(
"E field", always_available = .
true., updated_on_demand = .false., &
112 call this%quantities%add(
quantity_t(
"B field", always_available = .
true., updated_on_demand = .false., &
120 type(external_potential_t),
intent(inout) :: this
124 call this%deallocate_memory()
131 class(external_potential_t),
intent(inout) :: this
132 class(mesh_t),
intent(in) :: mesh
136 select case (this%type)
138 safe_allocate(this%pot(1:mesh%np))
140 safe_allocate(this%a_static(1:mesh%np, 1:this%space%dim))
142 if (this%space%periodic_dim < this%space%dim)
then
143 safe_allocate(this%pot(1:mesh%np))
144 safe_allocate(this%v_ext(1:mesh%np_part))
157 safe_deallocate_a(this%pot)
158 safe_deallocate_a(this%b_field)
159 safe_deallocate_a(this%a_static)
160 safe_deallocate_a(this%e_field)
161 safe_deallocate_a(this%v_ext)
173 select type (interaction)
177 message(1) =
"Unsupported interaction."
193 select type (interaction)
196 do ip = 1, interaction%system_np
197 interaction%partner_e_field(:, ip) = partner%e_field
198 interaction%partner_b_field(:, ip) =
m_zero
201 do ip = 1, interaction%system_np
202 interaction%partner_e_field(:, ip) =
m_zero
203 interaction%partner_b_field(:, ip) = partner%b_field
210 message(1) =
"Unsupported interaction."
222 class(
mesh_t),
intent(in) :: mesh
225 real(real64) :: pot_re, pot_im, r, xx(this%space%dim)
226 real(real64),
allocatable :: den(:), grx(:)
231 select case (this%type)
233 case (external_pot_usdef)
234 assert(
allocated(this%pot))
237 call mesh_r(mesh, ip, r, coords = xx)
239 this%pot(ip) = pot_re
243 assert(
allocated(this%pot))
247 write(
message(1),
'(a)')
'Error loading file '//trim(this%filename)//
'.'
248 write(
message(2),
'(a,i4)')
'Error code returned = ', err
253 assert(
allocated(this%pot))
255 safe_allocate(den(1:mesh%np))
258 call mesh_r(mesh, ip, r, coords = xx)
263 call dpoisson_solve(poisson, namespace, this%pot, den, all_nodes = .false.)
265 safe_deallocate_a(den)
268 assert(
allocated(this%b_field))
274 safe_allocate(grx(1:this%space%dim))
276 select case (this%space%dim)
278 select case (this%gauge_2d)
280 if (this%space%periodic_dim == 1)
then
281 message(1) =
"For 2D system, 1D-periodic, StaticMagneticField can only be "
282 message(2) =
"applied for StaticMagneticField2DGauge = linear_y."
287 grx(1:this%space%dim) = mesh%x(ip, 1:this%space%dim)
288 this%a_static(ip, 1) = -
m_half/
p_c * grx(2) * this%b_field(3)
289 this%a_static(ip, 2) =
m_half/
p_c * grx(1) * this%b_field(3)
294 grx(1:this%space%dim) = mesh%x(ip, 1:this%space%dim)
295 this%a_static(ip, 1) = -
m_one/
p_c * grx(2) * this%b_field(3)
296 this%a_static(ip, 2) =
m_zero
302 grx(1:this%space%dim) = mesh%x(ip, 1:this%space%dim)
303 this%a_static(ip, 1) = -
m_half/
p_c*(grx(2) * this%b_field(3) - grx(3) * this%b_field(2))
304 this%a_static(ip, 2) = -
m_half/
p_c*(grx(3) * this%b_field(1) - grx(1) * this%b_field(3))
305 this%a_static(ip, 3) = -
m_half/
p_c*(grx(1) * this%b_field(2) - grx(2) * this%b_field(1))
311 safe_deallocate_a(grx)
314 assert(
allocated(this%e_field))
316 if (this%space%periodic_dim < this%space%dim)
then
327 this%pot(ip) = sum(mesh%x(ip, this%space%periodic_dim + 1:this%space%dim) &
328 * this%e_field(this%space%periodic_dim + 1:this%space%dim))
332 this%v_ext(1:mesh%np) = this%pot(1:mesh%np)
333 do ip = mesh%np+1, mesh%np_part
334 this%v_ext(ip) = sum(mesh%x(ip, this%space%periodic_dim + 1:this%space%dim) &
335 * this%e_field(this%space%periodic_dim + 1:this%space%dim))
348 integer :: n_pot_block, row, read_data
352 integer :: dim, periodic_dim, idir
398 if (
parse_block(namespace,
'StaticExternalPotentials', blk) == 0)
then
401 do row = 0, n_pot_block-1
406 assert(read_data > 0)
408 call external_potentials%add(pot)
418 call parse_variable(namespace,
'PeriodicDimensions', 0, periodic_dim)
442 if (
parse_block(namespace,
'StaticMagneticField', blk) == 0)
then
461 call parse_variable(namespace,
'StaticMagneticField2DGauge', 0, pot%gauge_2d)
466 safe_allocate(pot%b_field(1:3))
474 if (periodic_dim == 2)
then
475 message(1) =
"StaticMagneticField cannot be applied in a 2D, 2D-periodic system."
478 if (pot%b_field(1)**2 + pot%b_field(2)**2 >
m_zero)
then
485 if (periodic_dim >= 2)
then
486 message(1) =
"In 3D, StaticMagneticField cannot be applied when the system is 2D- or 3D-periodic."
488 else if (periodic_dim == 1 .and. any(abs(pot%b_field(2:3)) >
m_zero))
then
489 message(1) =
"In 3D, 1D-periodic, StaticMagneticField must be zero in the y- and z-directions."
498 call external_potentials%add(pot)
516 if (
parse_block(namespace,
'StaticElectricField', blk) == 0)
then
522 safe_allocate(pot%e_field(1:dim))
527 if (idir <= periodic_dim .and. abs(pot%e_field(idir)) >
m_epsilon)
then
528 message(1) =
"Applying StaticElectricField in a periodic direction is only accurate for large supercells."
535 call external_potentials%add(pot)
548 type(
block_t),
intent(in) :: blk
549 integer,
intent(in) :: row
550 integer,
intent(out) :: read_data
552 integer :: ncols, icol, flag
566 if (pot%type >= 0)
then
567 message(1) =
'Error in reading the ExternalPotentials block'
577 call messages_input_error(namespace,
'ExternalPotentials',
"Unknown type of external potential")
584 if (icol >= ncols)
exit
590 case (option__staticexternalpotentials__file)
594 case (option__staticexternalpotentials__potential_formula)
599 if (pot%type /= external_pot_usdef)
then
600 call messages_input_error(namespace,
'ExternalPotentials',
'potential_formula can only be used with user_defined')
603 case (option__staticexternalpotentials__density_formula)
609 call messages_input_error(namespace,
'ExternalPotentials',
'density_formula can only be used with charge_density')
622 if (pot%type == external_pot_usdef .and. &
623 .not.
parameter_defined(option__staticexternalpotentials__potential_formula))
then
624 call messages_input_error(namespace,
'ExternalPotentials',
"The 'potential_formula' parameter is missing.")
629 call messages_input_error(namespace,
'ExternalPotentials',
"The 'density_formula' parameter is missing.")
643 integer(int64),
intent(in) :: param
657 integer(int64),
intent(in) :: param
662 call messages_input_error(namespace,
'ExternalPotentials',
"Duplicated parameter in external potential.")
subroutine check_duplication(param)
logical function parameter_defined(param)
integer, parameter, public external_pot_from_file
potential, defined in a file
subroutine, public load_external_potentials(external_potentials, namespace)
class(external_potential_t) function, pointer external_potential_init(namespace)
subroutine read_from_block(pot, namespace, blk, row, read_data)
integer, parameter, public external_pot_charge_density
user-defined function for charge density
subroutine external_potential_finalize(this)
subroutine external_potential_calculate(this, namespace, mesh, poisson)
subroutine external_potential_init_interaction_as_partner(partner, interaction)
subroutine external_potential_deallocate(this)
subroutine external_potential_copy_quantities_to_interaction(partner, interaction)
integer, parameter, public external_pot_static_efield
Static electric field.
integer, parameter, public external_pot_static_bfield
Static magnetic field.
subroutine external_potential_allocate(this, mesh)
real(real64), parameter, public m_zero
real(real64), parameter, public m_epsilon
real(real64), parameter, public m_half
real(real64), parameter, public p_c
Electron gyromagnetic ratio, see Phys. Rev. Lett. 130, 071801 (2023)
real(real64), parameter, public m_one
This module implements a simple hash table for non-negative integer keys and integer values.
subroutine, public iihash_end(h)
Free a hash table.
subroutine, public iihash_insert(h, key, val)
Insert a (key, val) pair into the hash table h.
integer function, public iihash_lookup(h, key, found)
Look up a value in the hash table h. If found is present, it indicates if key could be found in the t...
subroutine, public iihash_init(h)
Initialize a hash table h.
integer, parameter, public lorentz_force
This module defines classes and functions for interaction partners.
subroutine, public dio_function_input(filename, namespace, space, mesh, ff, ierr, map)
Reads a mesh function from file filename, and puts it into ff. If the map argument is passed,...
This module defines the meshes, which are used in Octopus.
pure subroutine, public mesh_r(mesh, ip, rr, origin, coords)
return the distance to the origin for a given grid point
subroutine, public messages_not_implemented(feature, namespace)
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)
subroutine, public messages_input_error(namespace, var, details, row, column)
integer function, public parse_block(namespace, name, blk, check_varinfo_)
subroutine, public dpoisson_solve(this, namespace, pot, rho, all_nodes, kernel, reset)
Calculates the Poisson equation. Given the density returns the corresponding potential.
This module defines the quantity_t class and the IDs for quantities, which can be exposed by a system...
subroutine, public conv_to_c_string(str)
converts to c string
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.
type(unit_system_t), public units_inp
the units systems for reading and writing
abstract class for general interaction partners
surrogate interaction class to avoid circular dependencies between modules.
This class implements the iteration counter used by the multisystem algorithms. As any iteration coun...
Lorenz force between a systems of particles and an electromagnetic field.
Describes mesh distribution to nodes.
Systems (system_t) can expose quantities that can be used to calculate interactions with other system...