1 /*
2 [auto_generated]
3 libs/numeric/odeint/examples/molecular_dynamics_cells.cpp
4
5 [begin_description]
6 Molecular dynamics example with cells.
7 [end_description]
8
9 Copyright 2009-2012 Karsten Ahnert
10 Copyright 2009-2012 Mario Mulansky
11
12 Distributed under the Boost Software License, Version 1.0.
13 (See accompanying file LICENSE_1_0.txt or
14 copy at http://www.boost.org/LICENSE_1_0.txt)
15 */
16
17 #include <boost/numeric/odeint.hpp>
18
19 #include <cstddef>
20 #include <vector>
21 #include <cmath>
22 #include <algorithm>
23 #include <tuple>
24 #include <iostream>
25 #include <random>
26
27 #include <boost/range/algorithm/for_each.hpp>
28 #include <boost/range/algorithm/sort.hpp>
29 #include <boost/range/algorithm/unique_copy.hpp>
30 #include <boost/range/algorithm_ext/iota.hpp>
31 #include <boost/iterator/zip_iterator.hpp>
32 #include <boost/iterator/transform_iterator.hpp>
33 #include <boost/iterator/permutation_iterator.hpp>
34 #include <boost/iterator/counting_iterator.hpp>
35
36 #include "point_type.hpp"
37
38
39
40
41
42
43
44
45 struct local_force
46 {
47 double m_gamma; // friction
local_forcelocal_force48 local_force( double gamma = 0.0 ) : m_gamma( gamma ) { }
49 template< typename Point >
operator ()local_force50 Point operator()( Point& x , Point& v ) const
51 {
52 return - m_gamma * v;
53 }
54 };
55
56
57 struct lennard_jones
58 {
59 double m_sigma;
60 double m_eps;
lennard_joneslennard_jones61 lennard_jones( double sigma = 1.0 , double eps = 0.1 ) : m_sigma( sigma ) , m_eps( eps ) { }
operator ()lennard_jones62 double operator()( double r ) const
63 {
64 double c = m_sigma / r;
65 double c3 = c * c * c;
66 double c6 = c3 * c3;
67 return 4.0 * m_eps * ( -12.0 * c6 * c6 / r + 6.0 * c6 / r );
68 }
69 };
70
71 template< typename F >
72 struct conservative_interaction
73 {
74 F m_f;
conservative_interactionconservative_interaction75 conservative_interaction( F const &f = F() ) : m_f( f ) { }
76 template< typename Point >
operator ()conservative_interaction77 Point operator()( Point const& x1 , Point const& x2 ) const
78 {
79 Point diff = x1 - x2;
80 double r = abs( diff );
81 double f = m_f( r );
82 return - diff / r * f;
83 }
84 };
85
86 template< typename F >
make_conservative_interaction(F const & f)87 conservative_interaction< F > make_conservative_interaction( F const &f )
88 {
89 return conservative_interaction< F >( f );
90 }
91
92
93
94
95
96
97
98 // force = interaction( x1 , x2 )
99 // force = local_force( x , v )
100 template< typename LocalForce , typename Interaction >
101 class md_system_bs
102 {
103 public:
104
105 typedef std::vector< double > vector_type;
106 typedef point< double , 2 > point_type;
107 typedef point< int , 2 > index_type;
108 typedef std::vector< point_type > point_vector;
109 typedef std::vector< index_type > index_vector;
110 typedef std::vector< size_t > hash_vector;
111 typedef LocalForce local_force_type;
112 typedef Interaction interaction_type;
113
114
115 struct params
116 {
117 size_t n;
118 size_t n_cell_x , n_cell_y , n_cells;
119 double x_max , y_max , cell_size;
120 double eps , sigma; // interaction strength, interaction radius
121 interaction_type interaction;
122 local_force_type local_force;
123 };
124
125
126 struct cell_functor
127 {
128 params const &m_p;
129
cell_functormd_system_bs::cell_functor130 cell_functor( params const& p ) : m_p( p ) { }
131
132 template< typename Tuple >
operator ()md_system_bs::cell_functor133 void operator()( Tuple const& t ) const
134 {
135 auto point = boost::get< 0 >( t );
136 size_t i1 = size_t( point[0] / m_p.cell_size ) , i2 = size_t( point[1] / m_p.cell_size );
137 boost::get< 1 >( t ) = index_type( i1 , i2 );
138 boost::get< 2 >( t ) = hash_func( boost::get< 1 >( t ) , m_p );
139 }
140 };
141
142
143
144 struct transform_functor
145 {
146 typedef size_t argument_type;
147 typedef size_t result_type;
148 hash_vector const* m_index;
transform_functormd_system_bs::transform_functor149 transform_functor( hash_vector const& index ) : m_index( &index ) { }
operator ()md_system_bs::transform_functor150 size_t operator()( size_t i ) const { return (*m_index)[i]; }
151 };
152
153
154
155 struct interaction_functor
156 {
157 hash_vector const &m_cells_begin;
158 hash_vector const &m_cells_end;
159 hash_vector const &m_order;
160 point_vector const &m_x;
161 point_vector const &m_v;
162 params const &m_p;
163 size_t m_ncellx , m_ncelly;
164
interaction_functormd_system_bs::interaction_functor165 interaction_functor( hash_vector const& cells_begin , hash_vector const& cells_end , hash_vector pos_order ,
166 point_vector const&x , point_vector const& v , params const &p )
167 : m_cells_begin( cells_begin ) , m_cells_end( cells_end ) , m_order( pos_order ) , m_x( x ) , m_v( v ) ,
168 m_p( p ) { }
169
170 template< typename Tuple >
operator ()md_system_bs::interaction_functor171 void operator()( Tuple const &t ) const
172 {
173 point_type x = periodic_bc( boost::get< 0 >( t ) , m_p ) , v = boost::get< 1 >( t );
174 index_type index = boost::get< 3 >( t );
175 size_t pos_hash = boost::get< 4 >( t );
176
177 point_type a = m_p.local_force( x , v );
178
179 for( int i=-1 ; i<=1 ; ++i )
180 {
181 for( int j=-1 ; j<=1 ; ++j )
182 {
183 index_type cell_index = index + index_type( i , j );
184 size_t cell_hash = hash_func( cell_index , m_p );
185 for( size_t ii = m_cells_begin[ cell_hash ] ; ii < m_cells_end[ cell_hash ] ; ++ii )
186 {
187 if( m_order[ ii ] == pos_hash ) continue;
188 point_type x2 = periodic_bc( m_x[ m_order[ ii ] ] , m_p );
189
190 if( cell_index[0] >= m_p.n_cell_x ) x2[0] += m_p.x_max;
191 if( cell_index[0] < 0 ) x2[0] -= m_p.x_max;
192 if( cell_index[1] >= m_p.n_cell_y ) x2[1] += m_p.y_max;
193 if( cell_index[1] < 0 ) x2[1] -= m_p.y_max;
194
195 a += m_p.interaction( x , x2 );
196 }
197 }
198 }
199 boost::get< 2 >( t ) = a;
200 }
201 };
202
203
204
205
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)206 md_system_bs( size_t n ,
207 local_force_type const& local_force = local_force_type() ,
208 interaction_type const& interaction = interaction_type() ,
209 double xmax = 100.0 , double ymax = 100.0 , double cell_size = 2.0 )
210 : m_p()
211 {
212 m_p.n = n;
213 m_p.x_max = xmax;
214 m_p.y_max = ymax;
215 m_p.interaction = interaction;
216 m_p.local_force = local_force;
217 m_p.n_cell_x = size_t( xmax / cell_size );
218 m_p.n_cell_y = size_t( ymax / cell_size );
219 m_p.n_cells = m_p.n_cell_x * m_p.n_cell_y;
220 m_p.cell_size = cell_size;
221 }
222
init_point_vector(point_vector & x) const223 void init_point_vector( point_vector &x ) const { x.resize( m_p.n ); }
224
operator ()(point_vector const & x,point_vector const & v,point_vector & a,double t) const225 void operator()( point_vector const& x , point_vector const& v , point_vector &a , double t ) const
226 {
227 // init
228 hash_vector pos_hash( m_p.n , 0 );
229 index_vector pos_index( m_p.n );
230 hash_vector pos_order( m_p.n , 0 );
231 hash_vector cells_begin( m_p.n_cells ) , cells_end( m_p.n_cells ) , cell_order( m_p.n_cells );
232
233 boost::iota( pos_order , 0 );
234 boost::iota( cell_order , 0 );
235
236 // calculate grid hash
237 // calcHash( m_dGridParticleHash, m_dGridParticleIndex, dPos, m_numParticles);
238 std::for_each(
239 boost::make_zip_iterator( boost::make_tuple( x.begin() , pos_index.begin() , pos_hash.begin() ) ) ,
240 boost::make_zip_iterator( boost::make_tuple( x.end() , pos_index.end() , pos_hash.end() ) ) ,
241 cell_functor( m_p ) );
242
243 // // sort particles based on hash
244 // // sortParticles(m_dGridParticleHash, m_dGridParticleIndex, m_numParticles);
245 boost::sort( pos_order , [&]( size_t i1 , size_t i2 ) -> bool {
246 return pos_hash[i1] < pos_hash[i2]; } );
247
248
249
250 // reorder particle arrays into sorted order and find start and end of each cell
251 std::for_each( cell_order.begin() , cell_order.end() , [&]( size_t i ) {
252 auto pos_begin = boost::make_transform_iterator( pos_order.begin() , transform_functor( pos_hash ) );
253 auto pos_end = boost::make_transform_iterator( pos_order.end() , transform_functor( pos_hash ) );
254 cells_begin[ i ] = std::distance( pos_begin , std::lower_bound( pos_begin , pos_end , i ) );
255 cells_end[ i ] = std::distance( pos_begin , std::upper_bound( pos_begin , pos_end , i ) );
256 } );
257
258 std::for_each(
259 boost::make_zip_iterator( boost::make_tuple(
260 x.begin() ,
261 v.begin() ,
262 a.begin() ,
263 pos_index.begin() ,
264 boost::counting_iterator< size_t >( 0 )
265 ) ) ,
266 boost::make_zip_iterator( boost::make_tuple(
267 x.end() ,
268 v.end() ,
269 a.end() ,
270 pos_index.end() ,
271 boost::counting_iterator< size_t >( m_p.n )
272 ) ) ,
273 interaction_functor( cells_begin , cells_end , pos_order , x , v , m_p ) );
274 }
275
bc(point_vector & x)276 void bc( point_vector &x )
277 {
278 for( size_t i=0 ; i<m_p.n ; ++i )
279 {
280 x[i][0] = periodic_bc( x[ i ][0] , m_p.x_max );
281 x[i][1] = periodic_bc( x[ i ][1] , m_p.y_max );
282 }
283 }
284
periodic_bc(double x,double xmax)285 static inline double periodic_bc( double x , double xmax )
286 {
287 double tmp = x - xmax * int( x / xmax );
288 return tmp >= 0.0 ? tmp : tmp + xmax;
289 }
290
291
periodic_bc(point_type const & x,params const & p)292 static inline point_type periodic_bc( point_type const& x , params const& p )
293 {
294 return point_type( periodic_bc( x[0] , p.x_max ) , periodic_bc( x[1] , p.y_max ) );
295 }
296
297
check_interval(int i,int max)298 static inline int check_interval( int i , int max )
299 {
300 int tmp = i % max;
301 return tmp >= 0 ? tmp : tmp + max;
302 }
303
304
hash_func(index_type index,params const & p)305 static inline size_t hash_func( index_type index , params const & p )
306 {
307 size_t i1 = check_interval( index[0] , p.n_cell_x );
308 size_t i2 = check_interval( index[1] , p.n_cell_y );
309 return i1 * p.n_cell_y + i2;
310 }
311
312 params m_p;
313 };
314
315
316 template< typename LocalForce , typename 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)317 md_system_bs< LocalForce , Interaction > make_md_system_bs( size_t n , LocalForce const &f1 , Interaction const &f2 ,
318 double xmax = 100.0 , double ymax = 100.0 , double cell_size = 2.0 )
319 {
320 return md_system_bs< LocalForce , Interaction >( n , f1 , f2 , xmax , ymax , cell_size );
321 }
322
323
324
325
326
327
328 using namespace boost::numeric::odeint;
329
330
331
main(int argc,char * argv[])332 int main( int argc , char *argv[] )
333 {
334 const size_t n1 = 32;
335 const size_t n2 = 32;
336 const size_t n = n1 * n2;
337 auto sys = make_md_system_bs( n , local_force() , make_conservative_interaction( lennard_jones() ) , 100.0 , 100.0 , 2.0 );
338 typedef decltype( sys ) system_type;
339 typedef system_type::point_vector point_vector;
340
341 std::mt19937 rng;
342 std::normal_distribution<> dist( 0.0 , 1.0 );
343
344 point_vector x , v;
345 sys.init_point_vector( x );
346 sys.init_point_vector( v );
347
348 for( size_t i=0 ; i<n1 ; ++i )
349 {
350 for( size_t j=0 ; j<n2 ; ++j )
351 {
352 size_t index = i * n2 + j;
353 x[index][0] = 10.0 + i * 2.0 ;
354 x[index][1] = 10.0 + j * 2.0 ;
355 v[index][0] = dist( rng ) ;
356 v[index][1] = dist( rng ) ;
357 }
358 }
359
360 velocity_verlet< point_vector > stepper;
361 const double dt = 0.025;
362 double t = 0.0;
363 // std::cout << "set term x11" << endl;
364 for( size_t oi=0 ; oi<10000 ; ++oi )
365 {
366 for( size_t ii=0 ; ii<50 ; ++ii,t+=dt )
367 stepper.do_step( sys , std::make_pair( std::ref( x ) , std::ref( v ) ) , t , dt );
368 sys.bc( x );
369
370 std::cout << "set size square" << "\n";
371 std::cout << "unset key" << "\n";
372 std::cout << "p [0:" << sys.m_p.x_max << "][0:" << sys.m_p.y_max << "] '-' pt 7 ps 0.5" << "\n";
373 for( size_t i=0 ; i<n ; ++i )
374 std::cout << x[i][0] << " " << x[i][1] << " " << v[i][0] << " " << v[i][1] << "\n";
375 std::cout << "e" << std::endl;
376 }
377
378 return 0;
379 }
380