Skip to content

Instantly share code, notes, and snippets.

@Warpten
Created May 15, 2017 20:25
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
Star You must be signed in to star a gist
Save Warpten/96e0bf6fd83a89222a7fcb801908868d to your computer and use it in GitHub Desktop.
#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;
}
}
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment