A [VMC](http://en.wikipedia.org/wiki/Variational_Monte_Carlo) implementation to find the ground state of the Simple Harmonic Oscillator.

sho.cc

// Compile with g++ -O2 -std=c++0x sho_test.cc -o sho_test
// or (gnu gcc 4.7 or above)
// Compile with g++ -O2 -std=c++11 sho_test.cc -o sho_test

// This package finds the ground state of the Simple Harmonic Oscillator
// with a test function \frac{\sqrt\alpha}{\pi^{\frac{1}{4}}} \exp( - \frac{1}{2} x^2\alpha^2) 

#include <cmath>
#include <iomanip>
#include <iostream>
#include <random>

namespace  helpers
{
    double probabiltyRatio(double x1, double x2, double alpha)
    {
        double prx1 = std::pow(std::exp( - 0.5 * x1 * x1 * alpha * alpha), 2);
        double prx2 = std::pow(std::exp( - 0.5 * x2 * x2 * alpha * alpha), 2);
        return (prx2 / prx1);
    }

    double localEnergy(double x, double alpha)
    {
        return alpha * alpha + x * x * (1 - std::pow(alpha, 4));
    }

    double stddev(double ex, double exsq)
    {
        return std::sqrt(exsq - std::pow(ex, 2));
    }
}


int main(int argc, char const *argv[])
{
    double alpha = 0.4;
    double mcs = 1000000;
    double step = 2.0;

    // Famous Merssene Twister engine, the best one :)
    // http://scicomp.stackexchange.com/questions/6886/hardware-random-number-generator-vs-pseudo-random-number-generator-in-the-battl
    std::mt19937_64 engine;

    // Useful distributions (in the non trivial way)
    // http://en.cppreference.com/w/cpp/numeric/random
    std::uniform_real_distribution<> zero_to_one(0.0,1.0);
    std::uniform_real_distribution<> minus_one_to_one(-1.0,1.0);

    while (alpha < 1.4)
    {
        double xNew;
        double x = minus_one_to_one(engine);
        double w;
        double s;
        double localEnergy;
        double energy = 0.0;
        double energysq = 0.0;

        int accepted = 0;

        for (int i = 0; i < mcs; ++i)
        {
            xNew = x + step * minus_one_to_one(engine);
            w = helpers::probabiltyRatio(x, xNew, alpha);
            s = zero_to_one(engine);
            if (w >= s) 
            {
                x = xNew;
                accepted++;
            }
            // You can print out the x values of some alphas so the probability
            // distribution of x ca be ploted in a histogram
            localEnergy = helpers::localEnergy(x, alpha);
            energy += localEnergy;
            energysq += localEnergy * localEnergy;
        }

        // Prints out: alpha, expected energy, and standard deviation of
        // the expected energy as well as the acceptance ration of the process
        // the "optimal" value of the aceptance ratio is about 0.5

        std::cout << alpha
                  << "    " << energy / mcs
                  << "    " << helpers::stddev(energy/mcs, energysq/mcs)
                  << "    " <<  (double) accepted / mcs << std::endl;
        alpha += 0.01;
    }
    return 0;
}