libs/numeric/odeint/examples/molecular_dynamics.cpp

   1 /*
   2  [auto_generated]
   3  libs/numeric/odeint/examples/molecular_dynamics.cpp
   4
   5  [begin_description]
   6  Molecular dynamics example.
   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 <vector>
  20 #include <iostream>
  21 #include <random>
  22
  23 using namespace boost::numeric::odeint;
  24
  25
  26
  27 using namespace std;
  28 #define tab "\t"
  29
  30 const size_t n1 = 16;
  31 const size_t n2 = 16;
  32
  33 struct md_system
  34 {
  35     static const size_t n = n1 * n2;
  36     typedef std::vector< double > vector_type;
  37
  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 )
  46     { }
  47
  48     static void init_vector_type( vector_type &x ) { x.resize( 2 * n ); }
  49
  50     void operator()( vector_type const& x , vector_type const& v , vector_type &a , double t ) const
  51     {
  52         for( size_t i=0 ; i<n ; ++i )
  53         {
  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 ] ;
  59         }
  60
  61         for( size_t i=0 ; i<n ; ++i )
  62         {
  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 )
  67             {
  68                 double xj = x[j] , yj = x[n+j];
  69                 xj = periodic_bc( xj , m_xmax );
  70                 yj = periodic_bc( yj , m_ymax );
  71
  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;
  79             }
  80         }
  81     }
  82
  83     void bc( vector_type &x )
  84     {
  85         for( size_t i=0 ; i<n ; ++i )
  86         {
  87             x[ i     ] = periodic_bc( x[ i     ] , m_xmax );
  88             x[ i + n ] = periodic_bc( x[ i + n ] , m_ymax );
  89         }
  90     }
  91
  92     inline double lennard_jones( double r ) const
  93     {
  94         double c = m_sigma / r;
  95         double c3 = c * c * c;
  96         double c6 = c3 * c3;
  97         return 4.0 * m_eps * ( -12.0 * c6 * c6 / r + 6.0 * c6 / r );
  98     }
  99
 100     static inline double periodic_bc( double x , double xmax )
 101     {
 102         return ( x < 0.0 ) ? x + xmax : ( x > xmax ) ? x - xmax : x ;
 103     }
 104
 105     double m_a;
 106     double m_gamma;
 107     double m_eps ;
 108     double m_sigma ;
 109     double m_xmax , m_ymax;
 110 };
 111
 112
 113
 114
 115
 116 int main( int argc , char *argv[] )
 117 {
 118     const size_t n = md_system::n;
 119     typedef md_system::vector_type vector_type;
 120
 121
 122     std::mt19937 rng;
 123     std::normal_distribution<> dist( 0.0 , 1.0 );
 124
 125     vector_type x , v;
 126     md_system::init_vector_type( x );
 127     md_system::init_vector_type( v );
 128
 129     for( size_t i=0 ; i<n1 ; ++i )
 130     {
 131         for( size_t j=0 ; j<n2 ; ++j )
 132         {
 133             x[i*n2+j  ] = 5.0 + i * 4.0 ;
 134             x[i*n2+j+n] = 5.0 + j * 4.0 ;
 135             v[i]   = dist( rng ) ;
 136             v[i+n] = dist( rng ) ;
 137         }
 138     }
 139
 140     velocity_verlet< vector_type > stepper;
 141     const double dt = 0.025;
 142     double t = 0.0;
 143     md_system sys;
 144     for( size_t oi=0 ; oi<100000 ; ++oi )
 145     {
 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 );
 148         sys.bc( x );
 149
 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;
 156     }
 157
 158
 159     return 0;
 160 }