1 /*
2  Copyright 2011 Mario Mulansky
3  Copyright 2012-2013 Karsten Ahnert
4 
5  Distributed under the Boost Software License, Version 1.0.
6  (See accompanying file LICENSE_1_0.txt or
7  copy at http://www.boost.org/LICENSE_1_0.txt)
8  */
9 
10 
11 /* strongly nonlinear hamiltonian lattice in 2d */
12 
13 #ifndef LATTICE2D_HPP
14 #define LATTICE2D_HPP
15 
16 #include <vector>
17 
18 #include <boost/math/special_functions/pow.hpp>
19 
20 using boost::math::pow;
21 
22 template< int Kappa , int Lambda >
23 struct lattice2d {
24 
25     const double m_beta;
26     std::vector< std::vector< double > > m_omega;
27 
lattice2dlattice2d28     lattice2d( const double beta )
29         : m_beta( beta )
30     { }
31 
32     template< class StateIn , class StateOut >
operator ()lattice2d33     void operator()( const StateIn &q , StateOut &dpdt )
34     {
35         // q and dpdt are 2d
36         const int N = q.size();
37 
38         int i;
39         for( i = 0 ; i < N ; ++i )
40         {
41             const int i_l = (i-1+N) % N;
42             const int i_r = (i+1) % N;
43             for( int j = 0 ; j < N ; ++j )
44             {
45             const int j_l = (j-1+N) % N;
46             const int j_r = (j+1) % N;
47             dpdt[i][j] = - m_omega[i][j] * pow<Kappa-1>( q[i][j] )
48                 - m_beta * pow<Lambda-1>( q[i][j] - q[i][j_l] )
49                 - m_beta * pow<Lambda-1>( q[i][j] - q[i][j_r] )
50                 - m_beta * pow<Lambda-1>( q[i][j] - q[i_l][j] )
51                 - m_beta * pow<Lambda-1>( q[i][j] - q[i_r][j] );
52             }
53         }
54     }
55 
56     template< class StateIn >
energylattice2d57     double energy( const StateIn &q , const StateIn &p )
58     {
59         // q and dpdt are 2d
60         const int N = q.size();
61         double energy = 0.0;
62         int i;
63         for( i = 0 ; i < N ; ++i )
64         {
65             const int i_l = (i-1+N) % N;
66             const int i_r = (i+1) % N;
67             for( int j = 0 ; j < N ; ++j )
68             {
69             const int j_l = (j-1+N) % N;
70             const int j_r = (j+1) % N;
71             energy += p[i][j]*p[i][j] / 2.0
72                         + m_omega[i][j] * pow<Kappa>( q[i][j] ) / Kappa
73                 + m_beta * pow<Lambda>( q[i][j] - q[i][j_l] ) / Lambda / 2
74                 + m_beta * pow<Lambda>( q[i][j] - q[i][j_r] ) / Lambda / 2
75                 + m_beta * pow<Lambda>( q[i][j] - q[i_l][j] ) / Lambda / 2
76                 + m_beta * pow<Lambda>( q[i][j] - q[i_r][j] ) / Lambda / 2;
77             }
78         }
79         return energy;
80     }
81 
82 
83     template< class StateIn , class StateOut >
local_energylattice2d84     double local_energy( const StateIn &q , const StateIn &p , StateOut &energy )
85     {
86         // q and dpdt are 2d
87         const int N = q.size();
88         double e = 0.0;
89         int i;
90         for( i = 0 ; i < N ; ++i )
91         {
92             const int i_l = (i-1+N) % N;
93             const int i_r = (i+1) % N;
94             for( int j = 0 ; j < N ; ++j )
95             {
96                 const int j_l = (j-1+N) % N;
97                 const int j_r = (j+1) % N;
98                 energy[i][j] = p[i][j]*p[i][j] / 2.0
99                     + m_omega[i][j] * pow<Kappa>( q[i][j] ) / Kappa
100                     + m_beta * pow<Lambda>( q[i][j] - q[i][j_l] ) / Lambda / 2
101                     + m_beta * pow<Lambda>( q[i][j] - q[i][j_r] ) / Lambda / 2
102                     + m_beta * pow<Lambda>( q[i][j] - q[i_l][j] ) / Lambda / 2
103                     + m_beta * pow<Lambda>( q[i][j] - q[i_r][j] ) / Lambda / 2;
104                 e += energy[i][j];
105             }
106         }
107         //rescale
108         e = 1.0/e;
109         for( i = 0 ; i < N ; ++i )
110             for( int j = 0 ; j < N ; ++j )
111                 energy[i][j] *= e;
112         return 1.0/e;
113     }
114 
load_potlattice2d115     void load_pot( const char* filename , const double W , const double gap ,
116                    const size_t dim )
117     {
118         std::ifstream in( filename , std::ios::in | std::ios::binary );
119         if( !in.is_open() ) {
120             std::cerr << "pot file not found: " << filename << std::endl;
121             exit(0);
122         } else {
123             std::cout << "using pot file: " << filename << std::endl;
124         }
125 
126         m_omega.resize( dim );
127         for( int i = 0 ; i < dim ; ++i )
128         {
129             m_omega[i].resize( dim );
130             for( size_t j = 0 ; j < dim ; ++j )
131             {
132                 if( !in.good() )
133                 {
134                     std::cerr << "I/O Error: " << filename << std::endl;
135                     exit(0);
136                 }
137                 double d;
138                 in.read( (char*) &d , sizeof(d) );
139                 if( (d < 0) || (d > 1.0) )
140                 {
141                     std::cerr << "ERROR: " << d << std::endl;
142                     exit(0);
143                 }
144                 m_omega[i][j] = W*d + gap;
145             }
146         }
147 
148     }
149 
generate_potlattice2d150     void generate_pot( const double W , const double gap , const size_t dim )
151     {
152         m_omega.resize( dim );
153         for( size_t i = 0 ; i < dim ; ++i )
154         {
155             m_omega[i].resize( dim );
156             for( size_t j = 0 ; j < dim ; ++j )
157             {
158                 m_omega[i][j] = W*static_cast<double>(rand())/RAND_MAX + gap;
159             }
160         }
161     }
162 
163 };
164 
165 #endif
166