]>
git.proxmox.com Git - ceph.git/blob - ceph/src/boost/libs/numeric/odeint/examples/molecular_dynamics.cpp
3 libs/numeric/odeint/examples/molecular_dynamics.cpp
6 Molecular dynamics example.
9 Copyright 2009-2012 Karsten Ahnert
10 Copyright 2009-2012 Mario Mulansky
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)
17 #include <boost/numeric/odeint.hpp>
23 using namespace boost::numeric::odeint
;
35 static const size_t n
= n1
* n2
;
36 typedef std::vector
< double > vector_type
;
38 md_system( double a
= 0.0 , // strength of harmonic oscillator
39 double gamma
= 0.0 , // friction
40 double eps
= 0.1 , // interaction strenght
41 double sigma
= 1.0 , // interaction radius
42 double xmax
= 150.0 , double ymax
= 150.0 )
43 : m_a( a
) , m_gamma( gamma
)
44 , m_eps( eps
) , m_sigma( sigma
)
45 , m_xmax( xmax
) , m_ymax( ymax
)
48 static void init_vector_type( vector_type
&x
) { x
.resize( 2 * n
); }
50 void operator()( vector_type
const& x
, vector_type
const& v
, vector_type
&a
, double t
) const
52 for( size_t i
=0 ; i
<n
; ++i
)
54 double diffx
= x
[i
] - 0.5 * m_xmax
, diffy
= x
[i
+n
] - 0.5 * m_ymax
;
55 double r2
= diffx
* diffx
+ diffy
* diffy
;
56 double r
= std::sqrt( r2
);
57 a
[ i
] = - m_a
* r
* diffx
- m_gamma
* v
[ i
] ;
58 a
[ n
+ i
] = - m_a
* r
* diffy
- m_gamma
* v
[ n
+ i
] ;
61 for( size_t i
=0 ; i
<n
; ++i
)
63 double xi
= x
[i
] , yi
= x
[n
+i
];
64 xi
= periodic_bc( xi
, m_xmax
);
65 yi
= periodic_bc( yi
, m_ymax
);
66 for( size_t j
=0 ; j
<i
; ++j
)
68 double xj
= x
[j
] , yj
= x
[n
+j
];
69 xj
= periodic_bc( xj
, m_xmax
);
70 yj
= periodic_bc( yj
, m_ymax
);
72 double diffx
= ( xj
- xi
) , diffy
= ( yj
- yi
);
73 double r
= sqrt( diffx
* diffx
+ diffy
* diffy
);
74 double f
= lennard_jones( r
);
75 a
[ i
] += diffx
/ r
* f
;
76 a
[ n
+ i
] += diffy
/ r
* f
;
77 a
[ j
] -= diffx
/ r
* f
;
78 a
[ n
+ j
] -= diffy
/ r
* f
;
83 void bc( vector_type
&x
)
85 for( size_t i
=0 ; i
<n
; ++i
)
87 x
[ i
] = periodic_bc( x
[ i
] , m_xmax
);
88 x
[ i
+ n
] = periodic_bc( x
[ i
+ n
] , m_ymax
);
92 inline double lennard_jones( double r
) const
94 double c
= m_sigma
/ r
;
95 double c3
= c
* c
* c
;
97 return 4.0 * m_eps
* ( -12.0 * c6
* c6
/ r
+ 6.0 * c6
/ r
);
100 static inline double periodic_bc( double x
, double xmax
)
102 return ( x
< 0.0 ) ? x
+ xmax
: ( x
> xmax
) ? x
- xmax
: x
;
109 double m_xmax
, m_ymax
;
116 int main( int argc
, char *argv
[] )
118 const size_t n
= md_system::n
;
119 typedef md_system::vector_type vector_type
;
123 std::normal_distribution
<> dist( 0.0 , 1.0 );
126 md_system::init_vector_type( x
);
127 md_system::init_vector_type( v
);
129 for( size_t i
=0 ; i
<n1
; ++i
)
131 for( size_t j
=0 ; j
<n2
; ++j
)
133 x
[i
*n2
+j
] = 5.0 + i
* 4.0 ;
134 x
[i
*n2
+j
+n
] = 5.0 + j
* 4.0 ;
136 v
[i
+n
] = dist( rng
) ;
140 velocity_verlet
< vector_type
> stepper
;
141 const double dt
= 0.025;
144 for( size_t oi
=0 ; oi
<100000 ; ++oi
)
146 for( size_t ii
=0 ; ii
<100 ; ++ii
,t
+=dt
)
147 stepper
.do_step( sys
, std::make_pair( std::ref( x
) , std::ref( v
) ) , t
, dt
);
150 std::cout
<< "set size square" << "\n";
151 std::cout
<< "unset key" << "\n";
152 std::cout
<< "p [0:" << sys
.m_xmax
<< "][0:" << sys
.m_ymax
<< "] '-' pt 7 ps 0.5" << "\n";
153 for( size_t i
=0 ; i
<n
; ++i
)
154 std::cout
<< x
[i
] << " " << x
[i
+n
] << " " << v
[i
] << " " << v
[i
+n
] << "\n";
155 std::cout
<< "e" << std::endl
;