変分モンテカルロ法のプログラム（C非標準のdrand48()を使ったC++版）

vmc_drand48.cpp

// 元のコード→ https://miyantarumi.hatenablog.com/entry/2022/04/08/080000
#include <array> // for std::array
#include <cmath> // for std::asin, std::cos, std::exp, std::fabs, std::hypot, std::log, std::sqrt, std::sin
#include <cstdint> // for std::int32_t
#include <cstdlib> // for drand48, srand48
#include <iomanip> // for std::setprecision
#include <ios> // for std::ios::fixed, std::ios::floatfield
#include <iostream> // for std::cout, std::endl
#include <utility> // for std::make_pair, std::pair
#include <vector> // for std::vector

namespace {
static auto constexpr ALPHA_0 = 0.8; // 変分を始める初期値
static auto constexpr H = 1.0E-3; // 数値微分の幅
static auto constexpr NSAMPLES = 30000;
static auto constexpr NSKIP = 4000;
static auto constexpr NWALKERS = 400;
static auto constexpr SEED = 20240102;
static auto constexpr SIGMA = 0.4; // metropoliswalkでの正規乱数の分散
static auto constexpr THRESHOLD = 1.0E-6;

class MyRand final {
public:
    MyRand(double mean, double sigma);
    MyRand() = delete;
    ~MyRand() = default;

    double normal_dist_rand();
    double rand() const;

private:
    double mean_;
    double sigma_;
};

struct Walker_pos {
    Walker_pos()
        : x_(NWALKERS, 0.0)
        , y_(NWALKERS, 0.0)
        , z_(NWALKERS, 0.0)
    {
    }
    std::vector<double> x_;
    std::vector<double> y_;
    std::vector<double> z_;
};

void initial_position(MyRand& mr, Walker_pos& wa_pos);
std::pair<double, double> metropolis_test(double alpha, MyRand& mr);
void run_vmc();
}

int main()
{
    run_vmc();

    return 0;
}

namespace {
inline MyRand::MyRand(double mean, double sigma)
    : mean_(mean)
    , sigma_(sigma)
{
    srand48(SEED);
}

inline double MyRand::normal_dist_rand()
{
    return sigma_ * std::sqrt(-2.0 * std::log(drand48()));
}

inline double MyRand::rand() const
{
    return drand48();
}

void initial_position(MyRand& mr, Walker_pos& wa_pos)
{
    for (auto i = 0; i < NWALKERS; i++) {
        wa_pos.x_[i] = (mr.rand() - 0.5) * 2.0 * 2.0;
        wa_pos.y_[i] = (mr.rand() - 0.5) * 2.0 * 2.0;
        wa_pos.z_[i] = (mr.rand() - 0.5) * 2.0 * 2.0;
    }
}

std::pair<double, double> metropolis_test(double alpha, MyRand& mr)
{
    std::vector<double> E_w(NWALKERS, 0.0);
    std::vector<double> Edlnpsi_w(NWALKERS, 0.0);
    std::vector<double> dlnpsi_w(NWALKERS, 0.0);

    Walker_pos wa_pos;

    initial_position(mr, wa_pos);

    for (auto i = 0; i < NSAMPLES; i++) {
        // Metropolis Testを行う

        for (auto j = 0; j < NWALKERS; j++) {
            auto const x = wa_pos.x_[j];
            auto const y = wa_pos.y_[j];
            auto const z = wa_pos.z_[j];

            auto const xp = x + mr.normal_dist_rand();
            auto const yp = y + mr.normal_dist_rand();
            auto const zp = z + mr.normal_dist_rand();         
            
            auto const r = std::hypot(x, y, z);
            auto const rp = std::hypot(xp, yp, zp);

            if (mr.rand() < std::exp(2.0 * alpha * (- rp + r))) {
                wa_pos.x_[j] = xp;
                wa_pos.y_[j] = yp;
                wa_pos.z_[j] = zp;
            }

            // Metropolis Testの終了
            // 変数を更新し終えた

            // walkerについて<E>や<Edlnpsi>,<dlnpsi>を計算する.
            if (i >= NSKIP) {
                auto const rtemp = std::hypot(wa_pos.x_[j], wa_pos.y_[j], wa_pos.z_[j]);
                E_w[j] += - 1.0 / rtemp - 0.5 * alpha * (alpha - 2.0 / rtemp);
                Edlnpsi_w[j] += -rtemp * (-1.0 / rtemp - 0.5 * alpha * (alpha - 2.0 / rtemp));
                dlnpsi_w[j] -= rtemp;
            }
        }
    }

    auto E = 0.0; // 変数の初期化
    auto Edlnpsi = 0.0;
    auto dlnpsi = 0.0;

    for (auto i = 0; i < NWALKERS; i++) {
        E_w[i] /= static_cast<double>(NSAMPLES - NSKIP);
        Edlnpsi_w[i] /= static_cast<double>(NSAMPLES - NSKIP);
        dlnpsi_w[i] /= static_cast<double>(NSAMPLES - NSKIP);
    }

    // walker全ての平均を取ることにする
    for (auto j = 0; j < NWALKERS; j++) {
        E += E_w[j];
        Edlnpsi += Edlnpsi_w[j];
        dlnpsi += dlnpsi_w[j];
    }

    E /= NWALKERS;
    Edlnpsi /= NWALKERS;
    dlnpsi /= NWALKERS; // dE/dalphaの計算をするための部品を計算した.

    auto const dE = 2.0 * (Edlnpsi - E * dlnpsi); // dE/dalphaの計算ができた.

    return std::make_pair(E, dE);
}

void run_vmc()
{
    MyRand mr(0.0, SIGMA);
    std::cout.setf(std::ios::fixed, std::ios::floatfield);

    double alpha_ini; // Newton法の始まり値

    auto alpha_fin = ALPHA_0; // Newton法の更新値

    do {
        alpha_ini = alpha_fin; // 前のループの結果を始まりの値に代入

        // **********************************************
        // alpha_iniでのエネルギーの微分dE/dalphaを計算する.
        // **********************************************

        auto const dE = metropolis_test(alpha_ini, mr).second;

        // ***************************
        // ddE/ddalphaの数値微分を行う.
        // **************************

        auto const alpha_plus = alpha_ini + H; // 微分のためにalphaを変化させる
        auto const alpha_minus = alpha_ini - H;

        // *********************************
        // *********************************
        // alpha_plusでのdE/dalphaを計算する.
        // *********************************
        // *********************************

        auto const dE_plus = metropolis_test(alpha_plus, mr).second;

        // *********************************
        // *********************************
        // alpha_minusでのdE/dalphaを計算する.
        // *********************************
        // *********************************

        auto const dE_minus = metropolis_test(alpha_minus, mr).second; // dE/dalpha_minusの計算ができた.

        // *************************************************
        // *************************************************
        // alphaをNewton法で求めるための部品をすべて計算し終えた
        // *************************************************
        // *************************************************

        // alphaの更新を行う

        auto const ddE = (dE_plus - dE_minus) / (2.0 * H); // dE/dalphaの数値微分ddEが計算できた.

        alpha_fin = alpha_ini - dE / ddE; // alphaを更新する.

        std::cout << std::setprecision(10)
                  << "alpha_fin = "
                  << alpha_fin
                  << ", alpha_ini = "
                  << alpha_ini
                  << '\n';

    } while (std::fabs(alpha_fin - alpha_ini) > THRESHOLD);

    auto const alpha_va = alpha_fin; // 変分で最適化した変分パラメータ

    std::cout << "エネルギーが最小になる変分パラメータはα = " << alpha_va << '\n';

    // ****************************************************
    // alpha_vaでの変分で求めた基底状態のエネルギーを計算しよう
    // ****************************************************

    std::cout << "基底状態のエネルギーは"
              << metropolis_test(alpha_va, mr).first
              << " (Hartree)"
              << std::endl;
}
}