/*
[auto_generated]
libs/numeric/odeint/examples/molecular_dynamics_cells.cpp
[begin_description]
Molecular dynamics example with cells.
[end_description]
Copyright 2009-2012 Karsten Ahnert
Copyright 2009-2012 Mario Mulansky
Distributed under the Boost Software License, Version 1.0.
(See accompanying file LICENSE_1_0.txt or
copy at http://www.boost.org/LICENSE_1_0.txt)
*/
#include <boost/numeric/odeint.hpp>
#include <cstddef>
#include <vector>
#include <cmath>
#include <algorithm>
#include <tuple>
#include <iostream>
#include <random>
#include <boost/range/algorithm/for_each.hpp>
#include <boost/range/algorithm/sort.hpp>
#include <boost/range/algorithm/unique_copy.hpp>
#include <boost/range/algorithm_ext/iota.hpp>
#include <boost/iterator/zip_iterator.hpp>
#include <boost/iterator/transform_iterator.hpp>
#include <boost/iterator/permutation_iterator.hpp>
#include <boost/iterator/counting_iterator.hpp>
#include "point_type.hpp"
struct local_force
{
double m_gamma; // friction
local_force( double gamma = 0.0 ) : m_gamma( gamma ) { }
template< typename Point >
Point operator()( Point& x , Point& v ) const
{
return - m_gamma * v;
}
};
struct lennard_jones
{
double m_sigma;
double m_eps;
lennard_jones( double sigma = 1.0 , double eps = 0.1 ) : m_sigma( sigma ) , m_eps( eps ) { }
double operator()( double r ) const
{
double c = m_sigma / r;
double c3 = c * c * c;
double c6 = c3 * c3;
return 4.0 * m_eps * ( -12.0 * c6 * c6 / r + 6.0 * c6 / r );
}
};
template< typename F >
struct conservative_interaction
{
F m_f;
conservative_interaction( F const &f = F() ) : m_f( f ) { }
template< typename Point >
Point operator()( Point const& x1 , Point const& x2 ) const
{
Point diff = x1 - x2;
double r = abs( diff );
double f = m_f( r );
return - diff / r * f;
}
};
template< typename F >
conservative_interaction< F > make_conservative_interaction( F const &f )
{
return conservative_interaction< F >( f );
}
// force = interaction( x1 , x2 )
// force = local_force( x , v )
template< typename LocalForce , typename Interaction >
class md_system_bs
{
public:
typedef std::vector< double > vector_type;
typedef point< double , 2 > point_type;
typedef point< int , 2 > index_type;
typedef std::vector< point_type > point_vector;
typedef std::vector< index_type > index_vector;
typedef std::vector< size_t > hash_vector;
typedef LocalForce local_force_type;
typedef Interaction interaction_type;
struct params
{
size_t n;
size_t n_cell_x , n_cell_y , n_cells;
double x_max , y_max , cell_size;
double eps , sigma; // interaction strength, interaction radius
interaction_type interaction;
local_force_type local_force;
};
struct cell_functor
{
params const &m_p;
cell_functor( params const& p ) : m_p( p ) { }
template< typename Tuple >
void operator()( Tuple const& t ) const
{
auto point = boost::get< 0 >( t );
size_t i1 = size_t( point[0] / m_p.cell_size ) , i2 = size_t( point[1] / m_p.cell_size );
boost::get< 1 >( t ) = index_type( i1 , i2 );
boost::get< 2 >( t ) = hash_func( boost::get< 1 >( t ) , m_p );
}
};
struct transform_functor : public std::unary_function< size_t , size_t >
{
hash_vector const* m_index;
transform_functor( hash_vector const& index ) : m_index( &index ) { }
size_t operator()( size_t i ) const { return (*m_index)[i]; }
};
struct interaction_functor
{
hash_vector const &m_cells_begin;
hash_vector const &m_cells_end;
hash_vector const &m_order;
point_vector const &m_x;
point_vector const &m_v;
params const &m_p;
size_t m_ncellx , m_ncelly;
interaction_functor( hash_vector const& cells_begin , hash_vector const& cells_end , hash_vector pos_order ,
point_vector const&x , point_vector const& v , params const &p )
: m_cells_begin( cells_begin ) , m_cells_end( cells_end ) , m_order( pos_order ) , m_x( x ) , m_v( v ) ,
m_p( p ) { }
template< typename Tuple >
void operator()( Tuple const &t ) const
{
point_type x = periodic_bc( boost::get< 0 >( t ) , m_p ) , v = boost::get< 1 >( t );
index_type index = boost::get< 3 >( t );
size_t pos_hash = boost::get< 4 >( t );
point_type a = m_p.local_force( x , v );
for( int i=-1 ; i<=1 ; ++i )
{
for( int j=-1 ; j<=1 ; ++j )
{
index_type cell_index = index + index_type( i , j );
size_t cell_hash = hash_func( cell_index , m_p );
for( size_t ii = m_cells_begin[ cell_hash ] ; ii < m_cells_end[ cell_hash ] ; ++ii )
{
if( m_order[ ii ] == pos_hash ) continue;
point_type x2 = periodic_bc( m_x[ m_order[ ii ] ] , m_p );
if( cell_index[0] >= m_p.n_cell_x ) x2[0] += m_p.x_max;
if( cell_index[0] < 0 ) x2[0] -= m_p.x_max;
if( cell_index[1] >= m_p.n_cell_y ) x2[1] += m_p.y_max;
if( cell_index[1] < 0 ) x2[1] -= m_p.y_max;
a += m_p.interaction( x , x2 );
}
}
}
boost::get< 2 >( t ) = a;
}
};
md_system_bs( size_t n ,
local_force_type const& local_force = local_force_type() ,
interaction_type const& interaction = interaction_type() ,
double xmax = 100.0 , double ymax = 100.0 , double cell_size = 2.0 )
: m_p()
{
m_p.n = n;
m_p.x_max = xmax;
m_p.y_max = ymax;
m_p.interaction = interaction;
m_p.local_force = local_force;
m_p.n_cell_x = size_t( xmax / cell_size );
m_p.n_cell_y = size_t( ymax / cell_size );
m_p.n_cells = m_p.n_cell_x * m_p.n_cell_y;
m_p.cell_size = cell_size;
}
void init_point_vector( point_vector &x ) const { x.resize( m_p.n ); }
void operator()( point_vector const& x , point_vector const& v , point_vector &a , double t ) const
{
// init
hash_vector pos_hash( m_p.n , 0 );
index_vector pos_index( m_p.n );
hash_vector pos_order( m_p.n , 0 );
hash_vector cells_begin( m_p.n_cells ) , cells_end( m_p.n_cells ) , cell_order( m_p.n_cells );
boost::iota( pos_order , 0 );
boost::iota( cell_order , 0 );
// calculate grid hash
// calcHash( m_dGridParticleHash, m_dGridParticleIndex, dPos, m_numParticles);
std::for_each(
boost::make_zip_iterator( boost::make_tuple( x.begin() , pos_index.begin() , pos_hash.begin() ) ) ,
boost::make_zip_iterator( boost::make_tuple( x.end() , pos_index.end() , pos_hash.end() ) ) ,
cell_functor( m_p ) );
// // sort particles based on hash
// // sortParticles(m_dGridParticleHash, m_dGridParticleIndex, m_numParticles);
boost::sort( pos_order , [&]( size_t i1 , size_t i2 ) -> bool {
return pos_hash[i1] < pos_hash[i2]; } );
// reorder particle arrays into sorted order and find start and end of each cell
std::for_each( cell_order.begin() , cell_order.end() , [&]( size_t i ) {
auto pos_begin = boost::make_transform_iterator( pos_order.begin() , transform_functor( pos_hash ) );
auto pos_end = boost::make_transform_iterator( pos_order.end() , transform_functor( pos_hash ) );
cells_begin[ i ] = std::distance( pos_begin , std::lower_bound( pos_begin , pos_end , i ) );
cells_end[ i ] = std::distance( pos_begin , std::upper_bound( pos_begin , pos_end , i ) );
} );
std::for_each(
boost::make_zip_iterator( boost::make_tuple(
x.begin() ,
v.begin() ,
a.begin() ,
pos_index.begin() ,
boost::counting_iterator< size_t >( 0 )
) ) ,
boost::make_zip_iterator( boost::make_tuple(
x.end() ,
v.end() ,
a.end() ,
pos_index.end() ,
boost::counting_iterator< size_t >( m_p.n )
) ) ,
interaction_functor( cells_begin , cells_end , pos_order , x , v , m_p ) );
}
void bc( point_vector &x )
{
for( size_t i=0 ; i<m_p.n ; ++i )
{
x[i][0] = periodic_bc( x[ i ][0] , m_p.x_max );
x[i][1] = periodic_bc( x[ i ][1] , m_p.y_max );
}
}
static inline double periodic_bc( double x , double xmax )
{
double tmp = x - xmax * int( x / xmax );
return tmp >= 0.0 ? tmp : tmp + xmax;
}
static inline point_type periodic_bc( point_type const& x , params const& p )
{
return point_type( periodic_bc( x[0] , p.x_max ) , periodic_bc( x[1] , p.y_max ) );
}
static inline int check_interval( int i , int max )
{
int tmp = i % max;
return tmp >= 0 ? tmp : tmp + max;
}
static inline size_t hash_func( index_type index , params const & p )
{
size_t i1 = check_interval( index[0] , p.n_cell_x );
size_t i2 = check_interval( index[1] , p.n_cell_y );
return i1 * p.n_cell_y + i2;
}
params m_p;
};
template< typename LocalForce , typename Interaction >
md_system_bs< LocalForce , Interaction > make_md_system_bs( size_t n , LocalForce const &f1 , Interaction const &f2 ,
double xmax = 100.0 , double ymax = 100.0 , double cell_size = 2.0 )
{
return md_system_bs< LocalForce , Interaction >( n , f1 , f2 , xmax , ymax , cell_size );
}
using namespace boost::numeric::odeint;
int main( int argc , char *argv[] )
{
const size_t n1 = 32;
const size_t n2 = 32;
const size_t n = n1 * n2;
auto sys = make_md_system_bs( n , local_force() , make_conservative_interaction( lennard_jones() ) , 100.0 , 100.0 , 2.0 );
typedef decltype( sys ) system_type;
typedef system_type::point_vector point_vector;
std::mt19937 rng;
std::normal_distribution<> dist( 0.0 , 1.0 );
point_vector x , v;
sys.init_point_vector( x );
sys.init_point_vector( v );
for( size_t i=0 ; i<n1 ; ++i )
{
for( size_t j=0 ; j<n2 ; ++j )
{
size_t index = i * n2 + j;
x[index][0] = 10.0 + i * 2.0 ;
x[index][1] = 10.0 + j * 2.0 ;
v[index][0] = dist( rng ) ;
v[index][1] = dist( rng ) ;
}
}
velocity_verlet< point_vector > stepper;
const double dt = 0.025;
double t = 0.0;
// std::cout << "set term x11" << endl;
for( size_t oi=0 ; oi<10000 ; ++oi )
{
for( size_t ii=0 ; ii<50 ; ++ii,t+=dt )
stepper.do_step( sys , std::make_pair( std::ref( x ) , std::ref( v ) ) , t , dt );
sys.bc( x );
std::cout << "set size square" << "\n";
std::cout << "unset key" << "\n";
std::cout << "p [0:" << sys.m_p.x_max << "][0:" << sys.m_p.y_max << "] '-' pt 7 ps 0.5" << "\n";
for( size_t i=0 ; i<n ; ++i )
std::cout << x[i][0] << " " << x[i][1] << " " << v[i][0] << " " << v[i][1] << "\n";
std::cout << "e" << std::endl;
}
return 0;
}
