boundstates1d.cpp

Go to the documentation of this file.

/**
 * @file boundstates1d.cpp
 * @brief Implements the boundstates1d physics/numerical component used by QPE experiments.
 * @ingroup qpe_quantum
 *
 * Source: @c core/quantum/boundstates1d.cpp. This file is part of Quantum Phenomena Explorer (QPE),
 * maintained by Officina Turini. See LICENSE_IT.txt / LICENSE.txt for the
 * QPE freeware source-available terms and THIRD_PARTY_NOTICES.txt for
 * third-party components.
 */
 
#include "boundstates1d.h"
 
#include "core/math/complexfield2d.h"
#include "core/quantum/potentialfield2d.h"
 
#include <algorithm>
#include <cmath>
#include <complex>
#include <cstring>
#include <limits>
#include <stdexcept>
 
namespace qpe {
namespace {
 
constexpr double kPi = 3.141592653589793238462643383279502884;
 
/**
 * @brief Performs finite positive.
 * @param value Value to assign or evaluate.
 */
bool finitePositive(double value) noexcept
{
    return std::isfinite(value) && value > 0.0;
}
 
/**
 * @brief Returns whether power of two.
 * @param value Value to assign or evaluate.
 */
bool isPowerOfTwo(std::size_t value) noexcept
{
    return value != 0u && (value & (value - 1u)) == 0u;
}
 
/**
 * @brief Performs fft1 d.
 * @param data The data value used by the operation.
 * @param inverse The inverse value used by the operation.
 */
void fft1D(std::vector<std::complex<double> > &data, bool inverse)
{
    const std::size_t count = data.size();
    if (!isPowerOfTwo(count)) {
        throw std::invalid_argument("BoundStates1D: X grid must be a power of two");
    }
 
    std::size_t j = 0u;
    for (std::size_t i = 1u; i < count; ++i) {
        std::size_t bit = count >> 1u;
        while (j & bit) {
            j ^= bit;
            bit >>= 1u;
        }
        j ^= bit;
        if (i < j) {
            std::swap(data[i], data[j]);
        }
    }
 
    const double sign = inverse ? 1.0 : -1.0;
    for (std::size_t length = 2u; length <= count; length <<= 1u) {
        const double angle = sign * 2.0 * kPi / static_cast<double>(length);
        const std::complex<double> wLength(std::cos(angle), std::sin(angle));
        const std::size_t half = length >> 1u;
        for (std::size_t block = 0u; block < count; block += length) {
            std::complex<double> w(1.0, 0.0);
            for (std::size_t offset = 0u; offset < half; ++offset) {
                const std::complex<double> even = data[block + offset];
                const std::complex<double> odd = data[block + offset + half] * w;
                data[block + offset] = even + odd;
                data[block + offset + half] = even - odd;
                w *= wLength;
            }
        }
    }
 
    if (inverse) {
        const double scale = 1.0 / static_cast<double>(count);
        for (std::complex<double> &value : data) {
            value *= scale;
        }
    }
}
 
/**
 * @brief Performs weighted norm.
 * @param v The v value used by the operation.
 * @param dx The dx value used by the operation.
 */
double weightedNorm(const std::vector<std::complex<double> > &v, double dx)
{
    double sum = 0.0;
    for (const std::complex<double> &value : v) {
        sum += std::norm(value);
    }
    return std::sqrt(std::max(0.0, sum * dx));
}
 
/**
 * @brief Returns normalize.
 * @param v The v value used by the operation.
 * @param dx The dx value used by the operation.
 */
bool normalize(std::vector<std::complex<double> > &v, double dx)
{
    const double n = weightedNorm(v, dx);
    if (!(n > 1.0e-30) || !std::isfinite(n)) {
        return false;
    }
    const double inv = 1.0 / n;
    for (std::complex<double> &value : v) {
        value *= inv;
    }
    return true;
}
 
std::complex<double> innerProduct(const std::vector<std::complex<double> > &a,
                                  const std::vector<std::complex<double> > &b,
                                  double dx)
{
    std::complex<double> sum(0.0, 0.0);
    for (std::size_t i = 0u; i < a.size(); ++i) {
        sum += std::conj(a[i]) * b[i];
    }
    return sum * dx;
}
 
/**
 * @brief Performs orthonormalize.
 * @param v The v value used by the operation.
 * @param basis The basis value used by the operation.
 * @param dx The dx value used by the operation.
 */
bool orthonormalize(std::vector<std::complex<double> > &v,
                    const std::vector<std::vector<std::complex<double> > > &basis,
                    double dx)
{
    // Two passes make the deflation noticeably more stable for higher modes.
    for (int pass = 0; pass < 2; ++pass) {
        for (const std::vector<std::complex<double> > &q : basis) {
            const std::complex<double> projection = innerProduct(q, v, dx);
            for (std::size_t i = 0u; i < v.size(); ++i) {
                v[i] -= projection * q[i];
            }
        }
    }
    return normalize(v, dx);
}
 
std::uint32_t nextRandom(std::uint32_t &state) noexcept
{
    // Deterministic xorshift32: enough for a reproducible non-symmetric trial
    // vector; this is not used for physical measurement sampling.
    state ^= state << 13;
    state ^= state >> 17;
    state ^= state << 5;
    return state;
}
 
std::vector<std::complex<double> > makeTrial(std::size_t count,
                                             const std::vector<double> &potential,
                                             double vMin,
                                             int stateIndex)
{
    std::vector<std::complex<double> > trial(count, std::complex<double>(0.0, 0.0));
    double spread = 0.0;
    for (double value : potential) {
        spread = std::max(spread, value - vMin);
    }
    const double scale = std::max(1.0, 0.25 * spread);
    std::uint32_t rng = 0x9E3779B9u ^ static_cast<std::uint32_t>(0x85EBCA6Bu * (stateIndex + 1));
 
    for (std::size_t i = 0u; i < count; ++i) {
        const double u = static_cast<double>(nextRandom(rng) & 0x00FFFFFFu) /
                         static_cast<double>(0x01000000u);
        const double randomValue = 2.0 * u - 1.0;
        const double weight = std::exp(-0.35 * (potential[i] - vMin) / scale);
        // Small deterministic harmonic component improves low-state convergence
        // while the random part prevents parity from hiding odd/even modes.
        const double harmonic = std::cos(2.0 * kPi * static_cast<double>(stateIndex + 1) *
                                         static_cast<double>(i) / static_cast<double>(count));
        trial[i] = weight * (0.75 * randomValue + 0.25 * harmonic);
    }
    return trial;
}
 
std::vector<double> spectralKinetic(std::size_t count,
                                    double dx,
                                    double mass,
                                    double hbar)
{
    std::vector<double> kinetic(count, 0.0);
    const double length = static_cast<double>(count) * dx;
    for (std::size_t k = 0u; k < count; ++k) {
        const long signedIndex = (k <= count / 2u)
            ? static_cast<long>(k)
            : static_cast<long>(k) - static_cast<long>(count);
        const double waveNumber = 2.0 * kPi * static_cast<double>(signedIndex) / length;
        kinetic[k] = hbar * hbar * waveNumber * waveNumber / (2.0 * mass);
    }
    return kinetic;
}
 
/**
 * @brief Performs expectation energy.
 * @param psi The psi value used by the operation.
 * @param potential Potential field used by the operation.
 * @param kinetic The kinetic value used by the operation.
 * @param dx The dx value used by the operation.
 */
double expectationEnergy(const std::vector<std::complex<double> > &psi,
                         const std::vector<double> &potential,
                         const std::vector<double> &kinetic,
                         double dx)
{
    std::vector<std::complex<double> > kPsi = psi;
    fft1D(kPsi, false);
    for (std::size_t i = 0u; i < kPsi.size(); ++i) {
        kPsi[i] *= kinetic[i];
    }
    fft1D(kPsi, true);
 
    std::complex<double> e(0.0, 0.0);
    for (std::size_t i = 0u; i < psi.size(); ++i) {
        const std::complex<double> hPsi = kPsi[i] + potential[i] * psi[i];
        e += std::conj(psi[i]) * hPsi;
    }
    return (e * dx).real();
}
 
/**
 * @brief Returns count nodes.
 * @param amplitude The amplitude value used by the operation.
 */
int countNodes(const std::vector<float> &amplitude)
{
    float maximum = 0.0f;
    for (float value : amplitude) {
        maximum = std::max(maximum, std::abs(value));
    }
    if (!(maximum > 0.0f)) {
        return 0;
    }
    const float threshold = maximum * 0.015f;
    int previousSign = 0;
    int nodes = 0;
    for (float value : amplitude) {
        if (std::abs(value) < threshold) {
            continue;
        }
        const int sign = value >= 0.0f ? 1 : -1;
        if (previousSign != 0 && sign != previousSign) {
            ++nodes;
        }
        previousSign = sign;
    }
    return nodes;
}
 
/**
 * @brief Performs edge mean.
 * @param values The values value used by the operation.
 * @param left The left value used by the operation.
 */
double edgeMean(const std::vector<double> &values, bool left)
{
    if (values.empty()) {
        return 0.0;
    }
    const std::size_t count = std::max<std::size_t>(4u, values.size() / 10u);
    double sum = 0.0;
    if (left) {
        for (std::size_t i = 0u; i < count; ++i) {
            sum += values[i];
        }
    } else {
        for (std::size_t i = values.size() - count; i < values.size(); ++i) {
            sum += values[i];
        }
    }
    return sum / static_cast<double>(count);
}
 
/**
 * @brief Performs finite vector.
 * @param values The values value used by the operation.
 */
bool finiteVector(const std::vector<double> &values) noexcept
{
    for (double value : values) {
        if (!std::isfinite(value)) {
            return false;
        }
    }
    return true;
}
 
} // namespace
 
bool BoundStates1D::isYInvariant(const PotentialField2D &potential,
                                 double tolerance) noexcept
{
    if (potential.empty() || potential.width() == 0u || potential.height() == 0u ||
        !std::isfinite(tolerance) || tolerance < 0.0) {
        return false;
    }
 
    const std::size_t referenceY = potential.height() / 2u;
    double scale = 1.0;
    for (std::size_t x = 0u; x < potential.width(); ++x) {
        scale = std::max(scale, std::abs(static_cast<double>(potential(x, referenceY))));
    }
    const double allowed = tolerance * scale;
 
    for (std::size_t y = 0u; y < potential.height(); ++y) {
        for (std::size_t x = 0u; x < potential.width(); ++x) {
            const double delta = std::abs(static_cast<double>(potential(x, y)) -
                                          static_cast<double>(potential(x, referenceY)));
            if (delta > allowed) {
                return false;
            }
        }
    }
    return true;
}
 
std::uint64_t BoundStates1D::potentialSignature(const PotentialField2D &potential) noexcept
{
    // FNV-1a over grid metadata and all float samples.  The signature is only
    // used to detect a changed preparation, never as a security primitive.
    std::uint64_t hash = 1469598103934665603ull;
    const auto mixBytes = [&hash](const void *data, std::size_t size) {
        const unsigned char *bytes = static_cast<const unsigned char *>(data);
        for (std::size_t i = 0u; i < size; ++i) {
            hash ^= static_cast<std::uint64_t>(bytes[i]);
            hash *= 1099511628211ull;
        }
    };
 
    const std::size_t width = potential.width();
    const std::size_t height = potential.height();
    const float dx = potential.dx();
    const float dy = potential.dy();
    mixBytes(&width, sizeof(width));
    mixBytes(&height, sizeof(height));
    mixBytes(&dx, sizeof(dx));
    mixBytes(&dy, sizeof(dy));
    if (!potential.empty()) {
        mixBytes(potential.data(), potential.size() * sizeof(PotentialField2D::Value));
    }
    return hash;
}
 
bool BoundStates1D::matches(const BoundStateSpectrum1D &spectrum,
                            const PotentialField2D &potential,
                            double mass,
                            double hbar) noexcept
{
    return spectrum.valid &&
           spectrum.sourceSignature == potentialSignature(potential) &&
           std::abs(spectrum.mass - mass) <= 1.0e-12 * std::max(1.0, std::abs(mass)) &&
           std::abs(spectrum.hbar - hbar) <= 1.0e-12 * std::max(1.0, std::abs(hbar)) &&
           spectrum.x.size() == potential.width();
}
 
BoundStateSpectrum1D BoundStates1D::solve(const PotentialField2D &potential,
                                          const BoundStateSolver1DParams &params)
{
    BoundStateSpectrum1D result;
    result.mass = params.mass;
    result.hbar = params.hbar;
    result.dx = potential.dx();
    result.sourceSignature = potentialSignature(potential);
 
    if (potential.empty() || potential.width() < 16u || potential.height() < 1u) {
        result.message = "Potential grid is empty or too small.";
        return result;
    }
    if (!isPowerOfTwo(potential.width())) {
        result.message = "Bound-state analysis requires a radix-2 X grid.";
        return result;
    }
    if (!finitePositive(params.mass) || !finitePositive(params.hbar) ||
        !finitePositive(params.imaginaryTimeStep) || params.maxIterations < 50 ||
        params.maxStates < 1 || params.maxStates > 16 ||
        !finitePositive(params.energyTolerance)) {
        result.message = "Invalid bound-state solver parameters.";
        return result;
    }
    if (!isYInvariant(potential, params.yInvarianceTolerance)) {
        result.message = "Potential is not invariant in Y. Use a well extended across Y for this 1D bound-state analysis.";
        return result;
    }
 
    const std::size_t n = potential.width();
    const std::size_t referenceY = potential.height() / 2u;
    result.x.resize(n);
    result.potential.resize(n);
    for (std::size_t i = 0u; i < n; ++i) {
        result.x[i] = potential.xCoordinate(i);
        result.potential[i] = static_cast<double>(potential(i, referenceY));
    }
    if (!finiteVector(result.potential)) {
        result.message = "Potential contains non-finite values.";
        return result;
    }
 
    const double dx = static_cast<double>(potential.dx());
    const double vMin = *std::min_element(result.potential.begin(), result.potential.end());
    result.continuumThreshold = std::min(edgeMean(result.potential, true),
                                         edgeMean(result.potential, false));
    const double potentialSpan = std::max(1.0,
        *std::max_element(result.potential.begin(), result.potential.end()) - vMin);
    const double boundMargin = 1.0e-5 * potentialSpan;
 
    const std::vector<double> kinetic = spectralKinetic(n, dx, params.mass, params.hbar);
    std::vector<double> halfPotentialDecay(n, 1.0);
    std::vector<double> kineticDecay(n, 1.0);
    for (std::size_t i = 0u; i < n; ++i) {
        halfPotentialDecay[i] = std::exp(-(result.potential[i] - vMin) *
                                         params.imaginaryTimeStep /
                                         (2.0 * params.hbar));
        kineticDecay[i] = std::exp(-kinetic[i] * params.imaginaryTimeStep / params.hbar);
    }
 
    std::vector<std::vector<std::complex<double> > > basis;
    basis.reserve(static_cast<std::size_t>(params.maxStates));
 
    for (int stateIndex = 0; stateIndex < params.maxStates; ++stateIndex) {
        std::vector<std::complex<double> > psi = makeTrial(n, result.potential, vMin, stateIndex);
        if (!orthonormalize(psi, basis, dx)) {
            break;
        }
 
        double previousEnergy = std::numeric_limits<double>::infinity();
        int stableChecks = 0;
        for (int iteration = 0; iteration < params.maxIterations; ++iteration) {
            for (std::size_t i = 0u; i < n; ++i) {
                psi[i] *= halfPotentialDecay[i];
            }
            fft1D(psi, false);
            for (std::size_t i = 0u; i < n; ++i) {
                psi[i] *= kineticDecay[i];
            }
            fft1D(psi, true);
            for (std::size_t i = 0u; i < n; ++i) {
                psi[i] *= halfPotentialDecay[i];
            }
            if (!orthonormalize(psi, basis, dx)) {
                break;
            }
 
            if ((iteration % 10) == 0) {
                const double energy = expectationEnergy(psi, result.potential, kinetic, dx);
                const double scale = std::max(1.0, std::abs(energy));
                if (std::isfinite(previousEnergy) &&
                    std::abs(energy - previousEnergy) <= params.energyTolerance * scale) {
                    ++stableChecks;
                } else {
                    stableChecks = 0;
                }
                previousEnergy = energy;
                if (iteration >= 120 && stableChecks >= 4) {
                    break;
                }
            }
        }
 
        if (!orthonormalize(psi, basis, dx)) {
            break;
        }
        const double energy = expectationEnergy(psi, result.potential, kinetic, dx);
        if (!std::isfinite(energy)) {
            break;
        }
 
        // Once the next orthogonal eigenmode reaches the edge continuum, stop:
        // higher modes are not spatially bound to the central well.
        if (!(energy < result.continuumThreshold - boundMargin)) {
            break;
        }
 
        // Rotate the global phase so the largest sample is real and positive.
        std::size_t peakIndex = 0u;
        double peakMagnitude = 0.0;
        for (std::size_t i = 0u; i < n; ++i) {
            const double magnitude = std::abs(psi[i]);
            if (magnitude > peakMagnitude) {
                peakMagnitude = magnitude;
                peakIndex = i;
            }
        }
        if (peakMagnitude > 0.0) {
            const double phase = std::arg(psi[peakIndex]);
            const std::complex<double> rotation(std::cos(-phase), std::sin(-phase));
            for (std::complex<double> &value : psi) {
                value *= rotation;
            }
        }
 
        BoundState1D state;
        state.quantumNumber = stateIndex + 1;
        state.energy = energy;
        state.amplitude.resize(n, 0.0f);
        for (std::size_t i = 0u; i < n; ++i) {
            state.amplitude[i] = static_cast<float>(psi[i].real());
        }
        state.nodeCount = countNodes(state.amplitude);
        result.states.push_back(state);
 
        // Rebuild the deflation basis from the displayed real representative.
        // For a real Hamiltonian this removes tiny numerical imaginary residue.
        std::vector<std::complex<double> > realState(n);
        for (std::size_t i = 0u; i < n; ++i) {
            realState[i] = static_cast<double>(state.amplitude[i]);
        }
        normalize(realState, dx);
        basis.push_back(realState);
    }
 
    result.valid = !result.states.empty();
    if (result.valid) {
        result.message = "Bound states computed.";
    } else {
        result.message = "No bound eigenstate was found below the edge continuum threshold.";
    }
    return result;
}
 
bool BoundStates1D::prepareEigenstate(ComplexField2D &field,
                                      const BoundStateSpectrum1D &spectrum,
                                      int stateIndex)
{
    if (!spectrum.valid || stateIndex < 0 ||
        stateIndex >= static_cast<int>(spectrum.states.size()) ||
        field.empty() || field.width() != spectrum.x.size()) {
        return false;
    }
 
    const BoundState1D &state = spectrum.states[static_cast<std::size_t>(stateIndex)];
    if (state.amplitude.size() != field.width()) {
        return false;
    }
 
    const double yMeasure = static_cast<double>(field.height()) * static_cast<double>(field.dy());
    if (!(yMeasure > 0.0)) {
        return false;
    }
    const float yScale = static_cast<float>(1.0 / std::sqrt(yMeasure));
    for (std::size_t y = 0u; y < field.height(); ++y) {
        for (std::size_t x = 0u; x < field.width(); ++x) {
            field(x, y) = ComplexField2D::Value(state.amplitude[x] * yScale, 0.0f);
        }
    }
    return field.normalize();
}
 
bool BoundStates1D::prepareSuperposition(ComplexField2D &field,
                                         const BoundStateSpectrum1D &spectrum,
                                         int stateIndexA,
                                         int stateIndexB,
                                         double relativePhaseRadians)
{
    if (!spectrum.valid || stateIndexA < 0 || stateIndexB < 0 ||
        stateIndexA == stateIndexB ||
        stateIndexA >= static_cast<int>(spectrum.states.size()) ||
        stateIndexB >= static_cast<int>(spectrum.states.size()) ||
        !std::isfinite(relativePhaseRadians) || field.empty() ||
        field.width() != spectrum.x.size()) {
        return false;
    }
 
    const BoundState1D &a = spectrum.states[static_cast<std::size_t>(stateIndexA)];
    const BoundState1D &b = spectrum.states[static_cast<std::size_t>(stateIndexB)];
    if (a.amplitude.size() != field.width() || b.amplitude.size() != field.width()) {
        return false;
    }
 
    const double yMeasure = static_cast<double>(field.height()) * static_cast<double>(field.dy());
    if (!(yMeasure > 0.0)) {
        return false;
    }
    const double oneOverSqrtTwo = 1.0 / std::sqrt(2.0);
    const std::complex<double> phase(std::cos(relativePhaseRadians),
                                     std::sin(relativePhaseRadians));
    const double yScale = 1.0 / std::sqrt(yMeasure);
 
    for (std::size_t y = 0u; y < field.height(); ++y) {
        for (std::size_t x = 0u; x < field.width(); ++x) {
            const std::complex<double> value = oneOverSqrtTwo *
                (static_cast<double>(a.amplitude[x]) +
                 phase * static_cast<double>(b.amplitude[x])) * yScale;
            field(x, y) = ComplexField2D::Value(static_cast<float>(value.real()),
                                                 static_cast<float>(value.imag()));
        }
    }
    return field.normalize();
}
 
} // namespace qpe