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