initial commit

.gitattributes

@@ -0,0 +1,2 @@
+* text=auto
+

.gitignore

@@ -0,0 +1,33 @@
+# Files generated by invoking Julia with --code-coverage
+*.jl.cov
+*.jl.*.cov
+
+# Files generated by invoking Julia with --track-allocation
+*.jl.mem
+
+# System-specific files and directories generated by the BinaryProvider and BinDeps packages
+# They contain absolute paths specific to the host computer, and so should not be committed
+deps/deps.jl
+deps/build.log
+deps/downloads/
+deps/usr/
+deps/src/
+
+# Build artifacts for creating documentation generated by the Documenter package
+docs/build/
+docs/site/
+
+# File generated by Pkg, the package manager, based on a corresponding Project.toml
+# It records a fixed state of all packages used by the project. As such, it should not be
+# committed for packages, but should be committed for applications that require a static
+# environment.
+Manifest.toml
+
+.env
+results/
+package-lock.json
+dist/
+node_modules/
+
+.vscode/
+

Project.toml

@@ -0,0 +1,9 @@
+name = "SchroedingerNumerics"
+uuid = "537d3920-04f8-45d5-b939-67e14d5b2f3b"
+authors = ["cnoelle <[email protected]>"]
+version = "0.1.0"
+
+[deps]
+JSON = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"

examples/harmonicOscillator/coherentStateClassical.json

@@ -0,0 +1,9 @@
+{
+    "type": "classical",
+    "scheme": {"id": "SymplecticEuler" , "mass": 1, "deltaX": 0.001},
+    "V_coefficients": [0,0,1],
+    "mass": 1.0,
+    "deltaT": 0.006283185307179587,
+    "t": 0.0,
+    "point": [1.0,0.0]
+}

src/ClassicalSystem.jl

@@ -0,0 +1,88 @@
+struct ClassicalSystem
+    initialState::AbstractPoint
+    hamiltonian::AbstractObservable
+    propagator::ClassicalPropagator
+    scheme::ClassicalNumericsScheme
+    deltaT::Real
+    # state
+    currentTime::Real
+    currentState::AbstractPoint
+    function ClassicalSystem(point::AbstractPoint, hamiltonian::AbstractObservable, deltaT::Real;
+            scheme::ClassicalNumericsScheme=SymplecticEuler(), t0::Real=0.)
+        propagator::ClassicalPropagator = incarnate(hamiltonian, scheme, deltaT)
+        return new(point, hamiltonian, propagator, scheme, deltaT, t0, point)
+    end # ClassicalSystem
+    function ClassicalSystem(other::ClassicalSystem, t::Real, point::AbstractPoint)
+        return new(other.initialState, other.hamiltonian, other.propagator, 
+            other.scheme, other.deltaT, t, point)
+    end # ClassicalSystem
+end # ClassicalSystem
+
+function scheme(system::ClassicalSystem)::NumericsScheme
+    return system.scheme
+end
+
+function hamiltonian(system::ClassicalSystem)::AbstractObservable
+    return system.hamiltonian
+end
+
+function deltaT(system::ClassicalSystem)::Real
+    return system.deltaT
+end
+
+function propagate(system::ClassicalSystem, timeSteps::Int=1)::ClassicalSystem
+    for _ in 1:timeSteps
+        point::AbstractPoint = propagateSingleTimestep(system.currentState, system.propagator)
+        system = ClassicalSystem(system, system.currentTime + system.deltaT, point)
+    end # for k
+    return system
+end # p
+
+function _writeSettings(system::ClassicalSystem, file::IO)
+    indent::String = "    "
+    write(file, "{\n")
+    write(file, indent, "\"type\": \"classical\",\n")
+    write(file, indent, "\"deltaT\": $(system.deltaT),\n")
+    write(file, indent, "\"scheme\": $(_schemeToJson(system.scheme))") 
+    try
+        _writePotentialJson(file, projectX(system.hamiltonian), 
+            indent=length(indent), startWithComma=true)
+    catch
+    end
+    write(file, "\n}\n")
+end # _writeSettings
+
+
+function trace(system::ClassicalSystem, timeSteps::Int=1000;
+        folder::String = "./results") # io::IO
+    Base.Filesystem.mkpath(folder)
+    settingsFile::String = joinpath(folder, "settings.json")
+    pointsFile::String = joinpath(folder, "points.csv")
+    open(settingsFile, "w") do settingsFile1
+        _writeSettings(system, settingsFile1)
+    end # settingsFile
+    V::Any = nothing
+    try
+        V = projectX(system.hamiltonian) # we expect this to be a function of x
+    catch
+    end
+    open(pointsFile, "w") do fileObservables
+        write(fileObservables, "x, p, E\n")
+        for _ in 1:timeSteps
+            # write expectation values
+            point::Point = pin(system.currentState)
+            xVal::Real = point.q
+            pVal::Real = point.p
+            energy::Real = system.hamiltonian(point)
+            write(fileObservables, "$(xVal), $(pVal), $(energy)\n")
+            system = propagate(system, 1)
+        end # for k
+    end # open fileObservables
+    return system
+
+end # trace
+
+export ClassicalSystem
+export propagate
+export trace
+

src/CrankNicolson.jl

@@ -0,0 +1,42 @@
+raw"Implements the Crank-Nicolson integration scheme for the Schrödinger equation, i.e.
+(1+i/(2\hbar) \hat H) \psi(t+1) = (1-i/(2\hbar) \hat H)\psi(t) "
+struct CrankNicolson <: QmNumericsScheme
+end # CrankNicolson
+
+function schemeId(::CrankNicolson)::String
+    return "CrankNicolson"
+end
+
+struct CNMatrix <: QuantumPropagator
+    representation::QmRepresentation
+    config::SchroedingerConfig
+    raw"(1+i*deltaT/(2\hbar) \hat H)"
+    A::AbstractArray{<:Complex, 2}
+    raw"(1-i*deltaT/(2\hbar) \hat H)"
+    B::AbstractArray{<:Complex, 2}
+    raw"The time evolution matrix (1+i/(2\hbar) \hat H)^(-1) * (1-i/(2\hbar) \hat H)"
+    AinvB::Union{AbstractArray{<:Complex, 2}, Nothing}
+end # CNMatrix
+
+"""
+Find the concrete matrix representation of the Hamiltonian for the selected grid
+"""
+function incarnate(hamiltonian::AbstractObservable, representation::QmRepresentation, 
+        scheme::CrankNicolson, deltaT::Real, invertMatrix::Bool=false)::QuantumPropagator
+    cfg = qmConfig(representation)
+    H = asOperator(hamiltonian, representation)
+    hbar = cfg.hbar
+    A = (LinearAlgebra.I + (im*deltaT/2/hbar) * H)
+    B = (LinearAlgebra.I - (im*deltaT/2/hbar) * H)
+    AinvB = invertMatrix ? LinearAlgebra.inv(A) * B : nothing
+    return CNMatrix(representation, cfg, A, B, AinvB)
+end # incarnate
+
+
+function propagateSingleTimestep(psi::AbstractWaveFunction, propagator::CNMatrix)::AbstractWaveFunction
+    valuesOld::AbstractArray{<:Complex, 1} = values(psi, propagator.representation)
+    valuesNew::AbstractArray{<:Complex, 1} = isnothing(propagator.AinvB) ? propagator.A\(propagator.B * valuesOld) : propagator.AinvB * valuesOld
+    return asWavefunction(valuesNew, propagator.representation)
+end # propagateSingleTimestep
+
+export CrankNicolson

src/DifferentiationUtils.jl

@@ -0,0 +1,102 @@
+# below: different types of operators associated to V, and helpers
+"
+Differentiation as a matrix operator
+
+* 1: forward/Newton (default)
+* 0: symmetric/midpoint
+* -1: backward
+"
+function _diffOperator(gridPoints::AbstractRange{<:Real}, order::Int; 
+        diffMethod::Int=1)::AbstractArray{<:Real, 2}
+    if order == 0
+        return LinearAlgebra.I
+    end # if
+    # In principle we can give a formula that works for all orders, see
+    # https://en.wikipedia.org/wiki/Numerical_differentiation#Higher_derivatives
+    # It is rather tricky, however, to take into account the boundary conditions, so we do it case by case 
+    numPoints::Int = length(gridPoints)
+    deltaX::Real = step(gridPoints)
+    if (order === 1)  # the first derivative may be calculated in different ways
+        # FIXME use ints?
+        if (diffMethod === 1)
+            diag1 = fill(-1., numPoints)
+            diag1[numPoints] = 0.
+            M1::AbstractArray{<:Real, 2} = LinearAlgebra.Bidiagonal(diag1, ones(numPoints-1), :U)
+            return M1/deltaX
+        elseif (diffMethod === 0)
+            upper = fill(1., numPoints-1)
+            lower = fill(-1., numPoints-1)
+            upper[1] = 2.
+            lower[numPoints-1] = -2.
+            diag2 = fill(0., numPoints)
+            diag2[1] = -2.
+            diag2[numPoints] = 2.
+            M2::AbstractArray{<:Real, 2} = LinearAlgebra.Tridiagonal(lower, diag2, upper)
+            return M2/deltaX/2
+        elseif (diffMethod === -1)
+            diag3 = fill(1., numPoints)
+            diag3[0] = 0.
+            M3::AbstractArray{<:Real, 2} = LinearAlgebra.Bidiagonal(-ones(numPoints), diag3, :L)
+            return M3/deltaX
+        end
+    elseif (order === 2)
+        diagonal::AbstractArray{Float64, 1} = fill(-2., numPoints)
+        diagonal[1] = -1.
+        diagonal[numPoints] = -1.
+        Msquared::AbstractArray{Float64, 2} = LinearAlgebra.Tridiagonal(ones(numPoints-1), diagonal, ones(numPoints-1))
+        return Msquared/(deltaX^2)
+    end # if
+    squares::Int = convert(Int, floor(order / 2))
+    M::AbstractArray{<:Real, 2} = _diffOperator(gridPoints, 2, diffMethod=diffMethod)^squares
+    return order % 2 === 0 ? M : M * _diffOperator(gridPoints, 1, diffMethod=diffMethod)
+end # _diffOperator
+
+"Differentiation as a matrix operator",
+function _diffOperator(gridPoints::AbstractArray{<:Real, 1}, order::Int; diffMethod::Int=1)::AbstractArray{<:Real, 2}
+    if order == 0
+        return LinearAlgebra.I
+    end # if
+    # this is the case that gridPoints is not equally spaced (is not a range)
+    # In principle we can give a formula that works for all orders, see
+    # https://en.wikipedia.org/wiki/Numerical_differentiation#Higher_derivatives
+    # It is rather tricky, however, to take into account the boundary conditions, so we do it case by case 
+    numPoints::Int = length(gridPoints)
+    if (order === 1)  # the first derivative may be calculated in different ways
+        if (diffMethod === 1)
+            upper1 = map(idx -> 1/(gridPoints[idx+1] - gridPoints[idx]), 1:(numPoints-1))
+            diag1 = -copy(upper1)
+            push!(diag1, 0.)
+            return LinearAlgebra.Bidiagonal(diag1, upper1, :U)
+        elseif (diffMethod === 0)
+            upper2 = map(idx -> 1/(gridPoints[idx+1] - gridPoints[idx-1]), 2:(numPoints-1))
+            lower2 = -copy(upper2)
+            diag2 = fill(0., numPoints-2)
+            firstUpper::Float64 = 1/(gridPoints[2] - gridPoints[1])
+            prepend!(upper2, firstUpper)
+            prepend!(diag2, -firstUpper)
+            lastLower::Float64 = 1/(gridPoints[numPoints] - gridPoints[numPoints-1])
+            push!(lower2, -lastLower)
+            push!(diag2, lastLower)
+            return LinearAlgebra.Tridiagonal(lower2, diag2, upper2)
+        elseif (diffMethod === -1)
+            diag3 = map(idx -> 1/(gridPoints[idx] - gridPoints[idx-1]), 2:numPoints)
+            lower3 = -copy(diag3)
+            prepend!(diag3, 0.)
+            return LinearAlgebra.Bidiagonal(diag3, lower3, :L)
+        end
+    elseif (order === 2)
+        upperSquare = map(idx -> 1/(gridPoints[idx+1] - gridPoints[idx])/(gridPoints[idx] - gridPoints[idx-1]), 2:(numPoints-1))
+        lowerSquare = copy(upperSquare)
+        diagSquare  = -2 * copy(upperSquare)
+        firstUpperSquare::Real = 1/(gridPoints[2] - gridPoints[1])^2
+        prepend!(upperSquare, firstUpperSquare)
+        prepend!(diagSquare, -firstUpperSquare)
+        lastLowerSquare::Real = 1/(gridPoints[numPoints] - gridPoints[numPoints-1])^2
+        push!(lowerSquare, lastLowerSquare)
+        push!(diagSquare, -lastLowerSquare)
+        return LinearAlgebra.Tridiagonal(lowerSquare, diagSquare, upperSquare)
+    end # if
+    squares::Int = convert(Int, floor(order / 2))
+    M::AbstractArray{<:Real, 2} = _diffOperator(gridPoints, 2, diffMethod=diffMethod)^squares
+    return order % 2 === 0 ? M : M * _diffOperator(gridPoints, 1, diffMethod=diffMethod)
+end # _diffOperator

src/ExactQuantumPropagation.jl

@@ -0,0 +1,31 @@
+"Get the next value of an exactly known wave function"
+struct ExactPropagation <: QmNumericsScheme
+end # ExactPropagation
+
+function schemeId(::ExactPropagation)::String
+    return "ExactQmPropagation"
+end
+
+struct ExactPropagator <: QuantumPropagator
+    deltaT::Real
+end # ExactPropagator
+
+function incarnate(hamiltonian::AbstractObservable, representation::QmRepresentation, 
+        scheme::ExactPropagation, deltaT::Real)::QuantumPropagator
+    return ExactPropagator(deltaT)
+end #specialize
+
+# propagate time step (qm)
+function propagateSingleTimestep(psi::ExactWaveFunction, propagator::ExactPropagator)::ExactWaveFunction
+    return ExactWaveFunction(psi.f, psi.timestamp + propagator.deltaT)
+end # propagateSingleTimestep
+
+function propagateSingleTimestep(psi::ExactSampledWaveFunction, propagator::ExactPropagator)::ExactSampledWaveFunction
+    return ExactSampledWaveFunction(psi.psi0, psi.f, psi.timestamp + propagator.deltaT)
+end # propagateSingleTimestep
+
+function propagateSingleTimestep(psi::TranslatedExactWaveFunction, propagator::ExactPropagator)::TranslatedExactWaveFunction
+    return TranslatedExactWaveFunction(psi.trajectory, propagateSingleTimestep(psi.psi, propagator), psi.hbar)
+end # propagateSingleTimestep
+
+export ExactPropagation

src/ExactTrajectory.jl

@@ -0,0 +1,41 @@
+struct ExactTrajectory <: ClassicalNumericsScheme
+    # A function Real -> Point (time -> phase space)
+    f::Any 
+end # struct
+
+function schemeId(scheme::ExactTrajectory)::String
+    return "ExactTrajectory"
+end
+
+struct ExactTrajectoryPropagator <: ClassicalPropagator
+    scheme::ExactTrajectory
+    deltaT::Real
+end # ExactTrajectoryPropagator
+
+function incarnate(hamiltonian::AbstractObservable, 
+        scheme::ExactTrajectory, deltaT::Real)::ExactTrajectoryPropagator
+    return ExactTrajectoryPropagator(scheme, deltaT)
+end #specialize
+
+struct TimedPoint <: AbstractPoint
+    q::Real
+    p::Real
+    t::Real
+end # TimedPoint
+
+function pin(point::TimedPoint)::Point
+    return Point(point.q, point.p)
+end # pin
+
+# propagate time step (classical)
+function propagateSingleTimestep(point::AbstractPoint, propagator::ExactTrajectoryPropagator)::TimedPoint
+    if !(point isa TimedPoint)
+        p = pin(point)
+        point = TimedPoint(p.q, p.p, 0.)
+    end # if
+    newT::Real = point.t + propagator.deltaT
+    newPoint::Point = propagator.scheme.f(newT)
+    return TimedPoint(newPoint.q, newPoint.p, newT)
+end # propagateSingleTimestep
+
+export ExactTrajectory