!> Author: Jabir Ali Ouassou
!> Category: Materials
!>
!> This module defines the data type 'material', which models the state of a physical material for a discretized range of
!> positions and energies. This is an abstract type, meaning that it is not intended to be instantiated on its own, but is
!> intended as a base type for physical materials like conductors, superconductors, and ferromagnets. In other words, this
!> type defines the essential data structures and program structure, while the derived subtypes will define actual physics.
module material_m
use :: stdio_m
use :: condmat_m
private
! Public declarations
public :: material_conf, material_load
! Type declarations
type, public, abstract :: material
! Material properties that affect the equilibrium state
real(wp) :: length = 1.00_wp !! Material length (L/ξ)
real(wp) :: thouless = 1.00_wp !! Thouless energy (ħD/L²)
real(wp) :: scattering = 0.01_wp !! Inelastic scattering (η/Δ₀)
logical :: transparent_a = .false. !! Interface transparency (left)
logical :: transparent_b = .false. !! Interface transparency (right)
logical :: phaselock = .false. !! Lock the center-of-mass phase?
! Material properties that affect the nonequilibrium state
logical :: nonequilibrium = .false. !! Equilibrium?
logical :: transverse = .false. !! Transverse potential gradients?
real(wp) :: voltage = 0.00_wp !! Voltage (eV/Δ₀)
real(wp) :: temperature = 0.01_wp !! Temperature (T/Tc)
real(wp) :: spinvoltage = 0.00_wp !! Spin-voltage (eVs/Δ₀)
real(wp) :: spintemperature = 0.00_wp !! Spin-temperature (Ts/Tc)
real(wp), dimension(1:3) :: spinaxis = [0,0,1] !! Spin quantization axis
! Physical state modelled using quasiclassical propagators
real(wp), allocatable :: energy(:) !! Energy domain
real(wp), allocatable :: location(:) !! Position domain
type(propagator), allocatable :: propagator(:,:) !! Propagator values
type(propagator), allocatable :: backup(:,:) !! Propagator backups
! Physical observables derived from the propagators above
real(wp), allocatable :: density(:,:,:) !! Spin-resolved density of states
real(wp), allocatable :: supercurrent(:,:) !! Charge, spin, heat, and spin-heat supercurrents
real(wp), allocatable :: lossycurrent(:,:) !! Charge, spin, heat, and spin-heat dissipative currents
real(wp), allocatable :: accumulation(:,:) !! Charge, spin, heat, and spin-heat accumulation
real(wp), allocatable :: magnetization(:,:) !! Magnetization due to exchange and Zeeman effects
complex(wp), allocatable :: correlation(:) !! Superconducting pair-correlations
! Hybrid structures are modeled as double-linked lists
integer :: order = 1 !! Simulation priority of this material
class(material), pointer :: material_a => null() !! Material to the left (default: vacuum)
class(material), pointer :: material_b => null() !! Material to the right (default: vacuum)
! Control parameters for the numerical solvers
integer :: iteration = 0 !! Used to keep track of selfconsistent iteration cycles
logical :: selfconsistent = .true. !! Whether to selfconsistently calculate the superconducting gap
logical :: boost = .true. !! Whether to use convergence acceleration methods
integer :: scaling = 128 !! Maximal mesh increase (range: 2^N, N>1)
integer :: method = 4 !! Runge—Kutta order (range: 2, 4, 6)
integer :: control = 2 !! Error control (1: defect, 2: global error, 3: 1 then 2, 4: 1 and 2)
real(wp) :: tolerance = 1e-10_wp !! Error tolerance (maximum defect or global error)
integer :: information = 0 !! Debug information (range: [-1,2])
real(wp) :: difference = 1e+10_wp !! Difference between iterations
! The following variables are used for input/output purposes, and should be modified by class(material) constructors
character(len=128) :: type_string = 'MATERIAL' !! Name of this material
contains
! These methods define how to update the physical state of the material
procedure(manipulate), deferred :: construct !! Construct object
procedure(initialize), deferred :: initialize !! Initialize object
procedure(manipulate), deferred :: update_prehook !! Executed before update
procedure(manipulate), deferred :: update_posthook !! Executed after update
procedure :: update => material_update !! Calculate propagators
procedure :: update_diffusion => diffusion_update !! Calculate propagators (in equilibrium)
procedure :: update_kinetic => kinetic_update !! Calculate propagators (nonequilibrium)
! These methods define the physical equations used by the update methods
procedure(diffusion_equation), deferred :: diffusion_equation !! Diffusion equation
procedure(diffusion_equation_a), deferred :: diffusion_equation_a !! Boundary condition (left)
procedure(diffusion_equation_b), deferred :: diffusion_equation_b !! Boundary condition (right)
procedure(kinetic_equation), deferred :: kinetic_equation !! Kinetic equation
procedure(kinetic_equation_a), deferred :: kinetic_equation_a !! Boundary condition (left)
procedure(kinetic_equation_b), deferred :: kinetic_equation_b !! Boundary condition (right)
! These methods define miscellaneous utility functions
procedure :: conf => material_conf !! Configures material parameters
procedure :: save => material_save !! Saves the state of the material
procedure :: load => material_load !! Loads the state of the material
end type
! Declare submodule procedures
interface
module subroutine kinetic_update(this)
class(material), target, intent(inout) :: this
end subroutine
module subroutine diffusion_update(this)
class(material), target, intent(inout) :: this
end subroutine
end interface
! Declare subclass procedures
abstract interface
subroutine initialize(this)
!! This interface is used for the deferred procedure initialize.
import material, wp
class(material), intent(inout) :: this
end subroutine
subroutine manipulate(this)
!! This interface is used for the deferred procedures construct, update_prehook, and update_posthook.
import material
class(material), intent(inout) :: this
end subroutine
pure subroutine diffusion_equation(this, p, e, z)
!! This interface is used for the deferred procedure diffusion_equation.
import material, propagator, wp
class(material), intent(in) :: this
type(propagator), intent(inout) :: p
complex(wp), intent(in) :: e
real(wp), intent(in) :: z
end subroutine
pure subroutine diffusion_equation_a(this, p, a, r, rt)
!! This interface is used for the deferred procedure diffusion_equation_a.
import material, propagator, spin, wp
class(material), intent(in) :: this
type(propagator), intent(in) :: p, a
type(spin), intent(inout) :: r, rt
end subroutine
pure subroutine diffusion_equation_b(this, p, b, r, rt)
!! This interface is used for the deferred procedure diffusion_equation_b.
import material, propagator, spin, wp
class(material), intent(in) :: this
type(propagator), intent(in) :: p, b
type(spin), intent(inout) :: r, rt
end subroutine
pure subroutine kinetic_equation(this, Gp, R, z)
!! This interface is used for the deferred procedure kinetic_equation.
import material, propagator, wp
class(material), intent(in) :: this
type(propagator), intent(in) :: Gp
complex(wp), dimension(0:7,0:7), intent(inout) :: R
real(wp), intent(in) :: z
end subroutine
pure subroutine kinetic_equation_a(this, Gp, Ga, Cp, Ca)
!! This interface is used for the deferred procedure kinetic_equation_a.
import material, propagator, wp
class(material), intent(in) :: this
type(propagator), intent(in) :: Gp, Ga
complex(wp), dimension(0:7,0:7), intent(out) :: Cp, Ca
end subroutine
pure subroutine kinetic_equation_b(this, Gp, Gb, Cp, Cb)
!! This interface is used for the deferred procedure kinetic_equation_b.
import material, propagator, wp
class(material), intent(in) :: this
type(propagator), intent(in) :: Gp, Gb
complex(wp), dimension(0:7,0:7), intent(out) :: Cp, Cb
end subroutine
end interface
contains
!--------------------------------------------------------------------------------!
! IMPLEMENTATION OF STATE UPDATE METHODS !
!--------------------------------------------------------------------------------!
subroutine material_update(this, bootstrap)
!! This subroutine updates the state of the material by solving the diffusion
!! equations for the equilibrium propagators, the kinetic equations for the
!! nonequilibrium propagators, and then calculating physical observables.
class(material), intent(inout) :: this !! Material that will be updated
logical, optional, intent(in) :: bootstrap !! Disable calculation of observables
! Only update materials with a positive order
if (this % order <= 0) then
return
end if
! Call the prehook method
call this % update_prehook
! Status information
if (this % information >= 0) then
block
character(len=10) :: str
write(str, '(2x,"[1m#",i0)') this % order
write(stdout,'(a,a,20x)') str, trim(this % type_string)
end block
end if
! Reset the difference since last update to zero
this % difference = 0.0_wp
! Solve the diffusion equations
call this % update_diffusion
! Solve the kinetic equations
if (this % nonequilibrium) then
call this % update_kinetic
end if
! Status information
if (this % information >= 0) then
write(stdout,'(6x,a,f10.8,a)') 'Max change: ',this % difference,' '
end if
! Stop here if bootstrapping
if (present(bootstrap)) then
if (bootstrap) then
return
end if
end if
! Call the posthook method
call this % update_posthook
end subroutine
!--------------------------------------------------------------------------------!
! IMPLEMENTATION OF UTILITY METHODS !
!--------------------------------------------------------------------------------!
subroutine material_conf(this, key, val)
!! Configure a material property based on a key-value pair.
use :: evaluate_m
class(material), intent(inout) :: this
character(*), intent(in) :: key
character(*), intent(in) :: val
select case(key)
case("length")
call evaluate(val, this%length)
this%thouless = 1/(this%length**2 + eps)
case("scattering")
call evaluate(val, this%scattering)
case("temperature")
call evaluate(val, this%temperature)
case("voltage")
call evaluate(val, this%voltage)
case("spinvoltage")
call evaluate(val, this%spinvoltage)
case("spintemperature")
call evaluate(val, this%spintemperature)
case("spinaxis")
call evaluate(val, this%spinaxis)
this%spinaxis = unitvector(this%spinaxis)
case("transverse")
call evaluate(val, this%transverse)
case("transparent_a")
call evaluate(val, this%transparent_a)
case("transparent_b")
call evaluate(val, this%transparent_b)
case("order")
call evaluate(val, this%order)
if (this%order > 16) then
call error("The order of the material must be be maximum 16.")
else if (this%order < 0) then
call error("The order of the material must be be minimum 0.")
end if
case("iteration")
call evaluate(val, this%iteration)
case("selfconsistent")
call evaluate(val, this%selfconsistent)
case("boost")
call evaluate(val, this%boost)
case("phaselock")
call evaluate(val, this%phaselock)
case("equilibrium")
call evaluate(val, this%nonequilibrium)
this % nonequilibrium = .not. this % nonequilibrium
case("nonequilibrium")
call evaluate(val, this%nonequilibrium)
case default
call error("Unknown material option '" // key // "'.")
end select
end subroutine
subroutine material_save(this)
!! Save a backup of the current material state.
class(material), intent(inout) :: this
! Make sure the backup exists
if (.not. allocated(this%backup)) then
allocate(this%backup(size(this%propagator,1),size(this%propagator,2)))
end if
! Make a backup of the propagators
call this % propagator % save(this % backup)
end subroutine
subroutine material_load(this)
!! Load a backup of a previous material state.
class(material), intent(inout) :: this
integer :: info
! Load saved propagators
call this % propagator % load(this % backup)
! Disable status messages
info = this%information
if (this%information >= 0) then
this%information = -1
end if
! Silently recalculate derived properties
call this % update_posthook
! Reenable status messages
this % information = info
end subroutine
end module
