gistfile1.txt

#include "box.hpp"
#include "object.hpp"
#include "predictor_corrector.hpp"
#include "calculator.hpp"
#include "extensions.hpp"
#include "logger.hpp"

#include <cassert>
#include <iostream>
#include <random>
#include <iomanip>

namespace simulation
{
    box::box() { }

    void box::energy(physics::energy::calculator* calculator)
    {
        _calculator = calculator;
    }

    void box::integrator(physics::trajectory::predictor_corrector* integrator)
    {
        _integrator = integrator;
    }

    physics::energy::calculator* box::energy()
    {
        return _calculator;
    }

    physics::trajectory::predictor_corrector* box::integrator()
    {
        return _integrator;
    }

    void box::set_density(f32 expected_density)
    {
        _density = expected_density;
    }

    void box::set_atom_mass(f32 mass)
    {
        particle_mass = mass;
    }

    void box::set_sigma(f32 sigma)
    {
        _density *= std::pow(sigma, 3.0f) / particle_mass;
    }

    void box::boundary_conditions()
    {
        for (object* obj : particles)
            obj->position -= std::round(obj->position / side_length) * side_length;
    }

    void box::nve(f32 time_step)
    {
        assert(integrator() != nullptr);
        assert(energy() != nullptr);

        integrator()->predict(time_step);

        energy()->calculate_potential(time_step);

        integrator()->correct(time_step);

        energy()->calculate_kinetic(time_step);
        energy()->correct();

        boundary_conditions();

        logging::energies << "Energy: " << energy()->kinetic_energy() << " " << energy()->potential_energy() << std::endl;
    }

    void box::print_positions() const
    {
        logging::trajectories << std::endl;
        // logging::trajectories << " BOX SIZE " << side_length << std::endl;
        // logging::trajectories << " BOUNDARIES [ -" << side_length / 2.0 << " " << side_length / 2.0 << " ]" << std::endl;
        // logging::trajectories << std::endl << " OBJECT POSITIONS" << std::endl;
        for (object* obj : particles)
            logging::trajectories << /*"      Ar   " <<*/ obj->position << std::endl;
        logging::trajectories << std::endl;
    }

    void box::generate_configuration(i32 count)
    {
        side_length = f32(std::pow(f32(count) / _density, 1.0 / 3.0));

        std::random_device r;
        std::mt19937 mt(r());
        std::uniform_real_distribution<> generator(-side_length / 2.0f, side_length / 2.0f);
        std::uniform_real_distribution<> varying_generator(1.01, 1.05);

        particles.resize(count);

        for (i32 i = 0; i < count; ++i)
            particles[i] = new simulation::object(generator(mt), generator(mt), generator(mt));

        for (i32 i = 0; i < count; ++i)
        {
            for (i32 j = 0; j < count; ++j)
            {
                if (i == j)
                    continue;

                math::vector<f64> vec = particles[j]->position - particles[i]->position;
                math::vector<f64> image_vec = vec - std::round(particles[j]->position / side_length) * side_length;
            
                f64 vec_length = image_vec.lengthSquared();
                if (vec_length < 1.0)
                {
                    f64 scaler = varying_generator(mt);
                    particles[j]->position += image_vec * scaler;
                    
                    j = -1;
                }
            }
        }

        for (i32 i = 0; i < count - 1; ++i)
            for (i32 j = i + 1; j < count; ++j)
                assert((particles[j]->position - particles[i]->position).length() < 1.0 && "Error in initial configuration!");

        boundary_conditions();
    }

    void box::print_pair_distances() const
    {
        logging::trajectories_extended << std::endl << " PAIR DISTANCES" << std::endl;
        logging::trajectories_extended << "          ";
        for (ui32 i = 0; i < particles.size(); ++i)
            logging::trajectories_extended << std::setw(10) << std::setfill(' ') << i;
        logging::trajectories_extended << std::endl;

        for (i32 i = i32(particles.size()) - 1; i >= 0; --i)
        {
            logging::trajectories_extended << std::setw(10) << std::setfill(' ') << i;
            for (i32 j = 0; j < i; ++j)
            {
                f64 dist = (particles[j]->position - particles[i]->position).length();
                logging::trajectories_extended << std::fixed << std::setprecision(6) << std::setw(10) << std::setfill(' ') << dist;
            }
            logging::trajectories_extended << std::endl;
        }

        logging::trajectories_extended << std::endl << " IMAGE PAIR DISTANCES" << std::endl;
        logging::trajectories_extended << "          ";
        for (ui32 i = 0; i < particles.size(); ++i)
            logging::trajectories_extended << std::setw(10) << std::setfill(' ') << i;
        logging::trajectories_extended << std::endl;

        for (i32 i = i32(particles.size()) - 1; i >= 0; --i)
        {
            logging::trajectories_extended << std::setw(10) << std::setfill(' ') << i;
            for (i32 j = 0; j < i; ++j)
            {
                math::vector<f64> dist = particles[j]->position - particles[i]->position;
                dist -= std::round(dist / side_length) * side_length;
                logging::trajectories_extended << std::fixed << std::setprecision(6) << std::setw(10) << std::setfill(' ') << dist.length();
            }
            logging::trajectories_extended << std::endl;
        }
    }
}