水素原子に対して、拡散モンテカルロ（DMC）法を用いて基底状態のエネルギーを求めるコード（Groq AI生成）

hydrogen_dmc.cpp

/*
    This code was generated by Groq AI and is released under the MIT License.

    MIT License

    Copyright (c) 2025 @dc1394

    Permission is hereby granted, free of charge, to any person obtaining a copy
    of this software and associated documentation files (the "Software"), to deal
    in the Software without restriction, including without limitation the rights
    to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
    copies of the Software, and to permit persons to whom the Software is
    furnished to do so, subject to the following conditions:

    The above copyright notice and this permission notice shall be included in all
    copies or substantial portions of the Software.

    THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
    IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
    FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
    AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
    LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
    OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
    SOFTWARE.
*/
#include <array>                        // for std::array
#include <boost/math/tools/minima.hpp>  // for boost::math::tools::brent_find_minima
#include <cmath>                        // for std::exp, std::log, std::sqrt
#include <cstdint>                      // for std::int32_t, std::size_t
#include <format>                       // for std::format
#include <iostream>                     // for std::cout, std::endl
#include <random>   // for std::mt19937, std::normal_distribution, std::random_device, std::uniform_real_distribution
#include <utility>  // for std::move, std::pair
#include <vector>   // for std::vector

namespace qmc {

// 3D位置を表す
using Position = std::array<double, 3>;

// ウォーカー構造
struct Walker
{
    Position pos;
    double   weight = 1.0;  // DMCで使用
};

// 試行波動関数関連
class TrialFunction
{
public:
    explicit TrialFunction(double alpha) : alpha_(alpha)
    {
    }

    // ψ(r) = exp(-α r)
    double psi(Position const & r) const
    {
        return std::exp(-alpha_ * norm(r));
    }

    // |ψ(r)|^2
    double psi_sq(Position const & r) const
    {
        auto const p = psi(r);
        return p * p;
    }

    // 局所エネルギー E_L(r) = -0.5 ∇²ψ/ψ + V(r)
    double local_energy(Position const & r) const
    {
        auto rr = norm(r);
        if (rr < 1e-10)
            rr = 1e-10;  // 特異点回避
        return -0.5 * alpha_ * alpha_ + (alpha_ - 1.0) / rr;
    }

    // ドリフト項 b = ∇ ln |ψ| = -α \vec{r} / r
    Position drift(Position const & r) const
    {
        auto const rr = norm(r);
        if (rr < 1e-10)
            return {0.0, 0.0, 0.0};
        auto const factor = -alpha_ / rr;  // For importance sampling, b = ∇ ln ψ
        return {factor * r[0], factor * r[1], factor * r[2]};
    }

private:
    double        alpha_;
    static double norm(const Position & r)
    {
        return std::sqrt(r[0] * r[0] + r[1] * r[1] + r[2] * r[2]);
    }
};

// VMCクラス
class VMC
{
public:
    VMC(double alpha, std::int32_t n_walkers, std::int32_t eq_steps, std::int32_t meas_steps, double delta)
        : tf_(alpha), n_walkers_(n_walkers), eq_steps_(eq_steps), meas_steps_(meas_steps), delta_(delta)
    {
        initialize_walkers();
    }

    // VMC実行、平均エネルギー + 標準誤差を返す
    std::pair<double, double> run()
    {
        equilibrate();
        return measure();
    }

    // 平衡化後のウォーカーを取得（DMC用）
    const std::vector<Walker> & get_equilibrated_walkers() const
    {
        return walkers_;
    }

private:
    TrialFunction       tf_;
    std::int32_t        n_walkers_, eq_steps_, meas_steps_;
    double              delta_;
    std::vector<Walker> walkers_;
    std::mt19937        rng_{std::random_device{}()};

    void initialize_walkers()
    {
        walkers_.resize(n_walkers_);
        std::normal_distribution<double> dist(0.0, 1.0);
        for (auto & w : walkers_) {
            w.pos = {dist(rng_), dist(rng_), dist(rng_)};
        }
    }

    void equilibrate()
    {
        for (auto step = 0; step < eq_steps_; ++step) {
            move_walkers();
        }
    }

    std::pair<double, double> measure()
    {
        std::vector<double> energies;
        energies.reserve(n_walkers_ * meas_steps_);
        auto total_moves    = 0;
        auto accepted_moves = 0;

        for (auto step = 0; step < meas_steps_; ++step) {
            accepted_moves += move_walkers();
            total_moves += n_walkers_;

            for (auto const & w : walkers_) {
                energies.push_back(tf_.local_energy(w.pos));
            }
        }
        auto const acceptance = static_cast<double>(accepted_moves) / total_moves;

        auto       sum = 0.0, sum_sq = 0.0;
        auto const n = energies.size();
        for (auto e : energies) {
            sum += e;
            sum_sq += e * e;
        }
        auto const mean   = sum / n;
        auto const var    = (sum_sq / n - mean * mean) / (n - 1.0);
        auto const stderr = std::sqrt(var / n);

        std::cout << std::format("VMC Acceptance rate: {:.7f}\n", acceptance);
        return {mean, stderr};
    }

    std::int32_t move_walkers()
    {
        std::uniform_real_distribution<double> uni(0.0, 1.0);
        std::normal_distribution<double>       norm(0.0, delta_);
        auto                                   accepted = 0;

        for (auto & w : walkers_) {
            auto new_pos = w.pos;
            for (auto d = 0; d < 3; ++d) {
                new_pos[d] += norm(rng_);
            }
            auto const ratio = tf_.psi_sq(new_pos) / tf_.psi_sq(w.pos);
            if (uni(rng_) < ratio) {
                w.pos = new_pos;
                ++accepted;
            }
        }
        return accepted;
    }
};

// DMCクラス
class DMC
{
public:
    DMC(double alpha, std::vector<Walker> && initial_walkers, double tau, std::int32_t n_steps, std::int32_t eq_steps,
        std::int32_t target_walkers)
        : tf_(alpha), tau_(tau), n_steps_(n_steps), eq_steps_(eq_steps), target_walkers_(target_walkers)
    {
        walkers_.reserve(target_walkers);
        for (auto const & iw : initial_walkers) {
            auto w   = iw;
            w.weight = 1.0;
            walkers_.push_back(w);
        }
        e_t_ = -0.5;  // 初期参照エネルギー
    }

    // DMC実行、平均エネルギー + 標準誤差を返す
    std::pair<double, double> run()
    {
        std::vector<double> energies;
        energies.reserve(n_steps_ - eq_steps_);
        std::vector<std::size_t> walker_counts;
        walker_counts.reserve(n_steps_);

        for (auto step = 0; step < n_steps_; ++step) {
            propagate();
            update_weights();
            branch();  // 今はweighted branching
            adjust_e_t();

            walker_counts.push_back(walkers_.size());

            if (step >= eq_steps_) {
                auto const avg_e = compute_average_energy();
                energies.push_back(avg_e);
            }

            if (step % 1000 == 0) {
                std::cout << std::format("DMC Step: {}, Walkers: {}, E_T: {:.7f}\n", step, walkers_.size(), e_t_);
            }
        }

        auto       sum = 0.0, sum_sq = 0.0;
        auto const n = energies.size();
        for (auto e : energies) {
            sum += e;
            sum_sq += e * e;
        }
        auto const mean   = sum / static_cast<double>(n);
        auto const var    = (sum_sq / static_cast<double>(n) - mean * mean) / (static_cast<double>(n) - 1.0);
        auto const stderr = std::sqrt(var / static_cast<double>(n));

        // ウォーカー数の統計
        auto       sum_w = 0.0, sum_w_sq = 0.0;
        auto const nw = walker_counts.size();
        for (auto wc : walker_counts) {
            sum_w += wc;
            sum_w_sq += wc * wc;
        }
        auto const mean_w = sum_w / static_cast<double>(nw);
        auto const var_w  = (sum_w_sq / static_cast<double>(nw) - mean_w * mean_w) / (static_cast<double>(nw) - 1.0);
        std::cout << std::format("DMC Average walkers: {:.7f} ± {:.7f} (target: {})\n", mean_w, std::sqrt(var_w),
                                 target_walkers_);

        return {mean, stderr};
    }

private:
    TrialFunction       tf_;
    double              tau_, e_t_;
    std::int32_t        n_steps_, eq_steps_, target_walkers_;
    std::vector<Walker> walkers_;
    std::mt19937        rng_{std::random_device{}()};

    void propagate()
    {
        std::normal_distribution<double> diff(0.0, std::sqrt(tau_));  // variance = τ (for D=1/2)

        for (auto & w : walkers_) {
            Position drift_move = tf_.drift(w.pos);
            for (std::int32_t d = 0; d < 3; ++d) {
                w.pos[d] += drift_move[d] * tau_ + diff(rng_);  // Δr = b τ + χ
            }
        }
    }

    void update_weights()
    {
        for (auto & w : walkers_) {
            auto const el = tf_.local_energy(w.pos);
            w.weight *= std::exp(-tau_ * (el - e_t_));
        }
    }

    void branch()
    {
        std::vector<Walker> new_walkers;
        new_walkers.reserve(walkers_.size() * 2);
        std::uniform_real_distribution<double> uni(0.0, 1.0);

        for (auto const & w : walkers_) {
            if (w.weight > 0.0) {
                auto const m = static_cast<std::int32_t>(w.weight + uni(rng_));
                for (auto i = 0; i < m; ++i) {
                    auto nw   = w;
                    nw.weight = 1.0;  // 重みリセット
                    new_walkers.push_back(nw);
                }
            }
        }
        walkers_ = std::move(new_walkers);

        // 死滅防止: 最小ウォーカー数
        if (walkers_.size() < 10) {
            std::cerr << "Warning: Walkers low, reinitializing..." << std::endl;
            // VMC-like reinitialization (simplified)
            std::normal_distribution<double> dist(0.0, 1.0);
            while (walkers_.size() < 10) {
                Walker w;
                w.pos    = {dist(rng_), dist(rng_), dist(rng_)};
                w.weight = 1.0;
                walkers_.push_back(w);
            }
        }
    }

    void adjust_e_t()
    {
        auto total_weight = 0.0;
        auto sum_el       = 0.0;
        for (auto const & w : walkers_) {
            total_weight += w.weight;
            sum_el += w.weight * tf_.local_energy(w.pos);
        }
        if (std::fabs(total_weight) < 1.0E-20) {
            return;
        }
        auto const avg_el = sum_el / total_weight;

        auto const log_term = std::log(static_cast<double>(target_walkers_) / walkers_.size());
        e_t_ *= 0.9 + (avg_el - 0.1 * log_term / tau_) * 0.1;  // スムーズ調整 (feedback constant 0.1)
    }

    double compute_average_energy() const
    {
        auto total_weight = 0.0;
        auto sum_el       = 0.0;
        for (auto const & w : walkers_) {
            total_weight += w.weight;
            sum_el += w.weight * tf_.local_energy(w.pos);
        }
        return (total_weight > 0.0) ? sum_el / total_weight : 0.0;
    }
};

}  // namespace qmc

int main()
{
    auto const n_walkers      = 1000;
    auto const eq_steps_vmc   = 1000;
    auto const meas_steps_vmc = 10000;
    auto const delta          = 1.0;

    auto const tau            = 0.001;  // 小さく調整
    auto const n_steps_dmc    = 10000;
    auto const eq_steps_dmc   = 1000;
    auto const target_walkers = 1000;

    // αの最適化 (VMCエネルギーを最小化)
    auto const vmc_energy_func = [&](double alpha) {
        qmc::VMC vmc(alpha, n_walkers, eq_steps_vmc, meas_steps_vmc, delta);
        auto [energy, _] = vmc.run();
        return energy;
    };

    auto const [optimal_alpha, min_energy] = boost::math::tools::brent_find_minima(vmc_energy_func, 0.5, 2.0, 32);

    std::cout << std::format("Optimal alpha: {:.7f}\n", optimal_alpha);

    // 最適αでVMC再実行（統計情報のため）
    qmc::VMC vmc_opt(optimal_alpha, n_walkers, eq_steps_vmc, meas_steps_vmc, delta);
    auto [vmc_energy, vmc_stderr] = vmc_opt.run();
    std::cout << std::format("VMC Energy: {:.7f} ± {:.7e}\n", vmc_energy, vmc_stderr);

    // DMC実行
    auto     init_walkers = vmc_opt.get_equilibrated_walkers();
    qmc::DMC dmc(optimal_alpha, std::move(init_walkers), tau, n_steps_dmc, eq_steps_dmc, target_walkers);
    auto const [dmc_energy, dmc_stderr] = dmc.run();
    std::cout << std::format("DMC Energy: {:.7f} ± {:.7e}\n", dmc_energy, dmc_stderr);

    return 0;
}

水素原子に対して、拡散モンテカルロ（DMC）法を用いて基底状態のエネルギーを求めるコード（Groq AI生成）です。
試行関数をΨ(r) = exp(-αr)とし、まず変分モンテカルロ（VMC）法でαを最適化してからDMCを行います。
コンパイルには、C++20に対応したコンパイラとBoost C++ライブラリが必要です。