1 /*
2 * gauss_packet.cpp
3 *
4 * Schroedinger equation with potential barrier and periodic boundary conditions
5 * Initial Gauss packet moving to the right
6 *
7 * pipe output into gnuplot to see animation
8 *
9 * Implementation of Hamilton operator via MTL library
10 *
11 * Copyright 2011-2013 Mario Mulansky
12 * Copyright 2011-2012 Karsten Ahnert
13 *
14 * Distributed under the Boost Software License, Version 1.0.
15 * (See accompanying file LICENSE_1_0.txt or
16 * copy at http://www.boost.org/LICENSE_1_0.txt)
17 */
18
19
20 #include <iostream>
21 #include <complex>
22
23 #include <boost/numeric/odeint.hpp>
24 #include <boost/numeric/odeint/external/mtl4/mtl4.hpp>
25
26 #include <boost/numeric/mtl/mtl.hpp>
27
28
29 using namespace std;
30 using namespace boost::numeric::odeint;
31
32 typedef mtl::dense_vector< complex< double > > state_type;
33
34 struct hamiltonian {
35
36 typedef mtl::compressed2D< complex< double > > matrix_type;
37 matrix_type m_H;
38
hamiltonianhamiltonian39 hamiltonian( const int N ) : m_H( N , N )
40 {
41 // constructor with zero potential
42 m_H = 0.0;
43 initialize_kinetic_term();
44 }
45
46 //template< mtl::compressed2D< double > >
hamiltonianhamiltonian47 hamiltonian( mtl::compressed2D< double > &V ) : m_H( num_rows( V ) , num_cols( V ) )
48 {
49 // use potential V in hamiltonian
50 m_H = complex<double>( 0.0 , -1.0 ) * V;
51 initialize_kinetic_term();
52 }
53
initialize_kinetic_termhamiltonian54 void initialize_kinetic_term( )
55 {
56 const int N = num_rows( m_H );
57 mtl::matrix::inserter< matrix_type , mtl::update_plus< complex<double> > > ins( m_H );
58 const double z = 1.0;
59 // fill diagonal and upper and lower diagonal
60 for( int i = 0 ; i<N ; ++i )
61 {
62 ins[ i ][ (i+1) % N ] << complex< double >( 0.0 , -z );
63 ins[ i ][ i ] << complex< double >( 0.0 , z );
64 ins[ (i+1) % N ][ i ] << complex< double >( 0.0 , -z );
65 }
66 }
67
operator ()hamiltonian68 void operator()( const state_type &psi , state_type &dpsidt , const double t )
69 {
70 dpsidt = m_H * psi;
71 }
72
73 };
74
75 struct write_for_gnuplot
76 {
77 size_t m_every , m_count;
78
write_for_gnuplotwrite_for_gnuplot79 write_for_gnuplot( size_t every = 10 )
80 : m_every( every ) , m_count( 0 ) { }
81
operator ()write_for_gnuplot82 void operator()( const state_type &x , double t )
83 {
84 if( ( m_count % m_every ) == 0 )
85 {
86 //clog << t << endl;
87 cout << "p [0:" << mtl::size(x) << "][0:0.02] '-'" << endl;
88 for( size_t i=0 ; i<mtl::size(x) ; ++i )
89 {
90 cout << i << "\t" << norm(x[i]) << "\n";
91 }
92 cout << "e" << endl;
93 }
94
95 ++m_count;
96 }
97 };
98
99 static const int N = 1024;
100 static const int N0 = 256;
101 static const double sigma0 = 20;
102 static const double k0 = -1.0;
103
main(int argc,char ** argv)104 int main( int argc , char** argv )
105 {
106 state_type x( N , 0.0 );
107
108 // initialize gauss packet with nonzero velocity
109 for( int i=0 ; i<N ; ++i )
110 {
111 x[i] = exp( -(i-N0)*(i-N0) / ( 4.0*sigma0*sigma0 ) ) * exp( complex< double >( 0.0 , k0*i ) );
112 //x[i] += 2.0*exp( -(i+N0-N)*(i+N0-N) / ( 4.0*sigma0*sigma0 ) ) * exp( complex< double >( 0.0 , -k0*i ) );
113 }
114 x /= mtl::two_norm( x );
115
116 typedef runge_kutta4< state_type > stepper;
117
118 // create potential barrier
119 mtl::compressed2D< double > V( N , N );
120 V = 0.0;
121 {
122 mtl::matrix::inserter< mtl::compressed2D< double > > ins( V );
123 for( int i=0 ; i<N ; ++i )
124 {
125 //ins[i][i] << 1E-4*(i-N/2)*(i-N/2);
126
127 if( i < N/2 )
128 ins[ i ][ i ] << 0.0 ;
129 else
130 ins[ i ][ i ] << 1.0 ;
131
132 }
133 }
134
135 // perform integration, output can be piped to gnuplot
136 integrate_const( stepper() , hamiltonian( V ) , x , 0.0 , 1000.0 , 0.1 , write_for_gnuplot( 10 ) );
137
138 clog << "Norm: " << mtl::two_norm( x ) << endl;
139
140 return 0;
141 }
142