Changed all literalincludes --> runblock / triqs_example

This is to avoid keeping code snippets that do not work in the doc. At least there will be an error message.

doc/reference/c++/arrays/h5_stack.rst

@@ -0,0 +1,51 @@
+.. highlight:: c
+
+h5::array_stack : stacking arrays or scalars 
+================================================================
+
+h5::array_stack writes a sequences of arrays of the same shape (or of scalars) into an hdf5 array with one more dimension, unlimited in the stacking direction.
+
+It is typically used to store a Monte-Carlo data series for later analysis.
+
+* If the base of the stack is an array of rank R, the resulting hdf5 array will be of rank R+1.
+
+* If the base of the stack is a simple number (double, int, ...), R=0.
+
+* The syntax is simple : 
+
+ * The << operator piles up an array/scalar onto the stack.
+ * The ++ operator advances by one slice in the stack.
+ * The () operator returns a view on the current slice of the stack.
+
+* The stack is bufferized in memory (`bufsize` parameter), so that the file access does not happen too often.
+
+* NB : beware to complex numbers ---> REF TO COMPLEX
+
+Reference 
+------------
+
+Here is the :doxy:`full C++ documentation<triqs::arrays::array_stack>` for this class.
+
+.. :
+    Breathe Documentation 
+  --------------------------
+  
+  .. doxygenclass:: triqs::arrays::array_stack
+   :project: arrays
+   :members:
+ 
+ 
+
+Tutorial
+-----------
+
+A simple example with a stack of double: 
+
+.. triqs_example:: examples_code/h5_stack_ex_sca.cpp
+
+A simple example with a stack of array of rank 2 :
+
+.. triqs_example:: examples_code/h5_stack_ex.cpp
+
+
+

doc/reference/c++/clef/examples/lazy_sum.rst

@@ -6,10 +6,4 @@ A lazy sum
 Here is a little functional `sum` that sums a function f over various domains
 and accepts lazy expressions as arguments.
 
-.. literalinclude:: src/sum_functional.cpp
-
-Compiling and running this code reads : 
-  
-.. literalinclude:: src/sum_functional.output
- 
- 
+.. triqs_example:: src/sum_functional.cpp

doc/reference/c++/statistics/ising2d.rst

@@ -4,8 +4,5 @@ Full example: Monte-Carlo simulation of the 2D Ising model
 ===========================================================
 
 
-.. literalinclude:: src/ising2d.cpp
+.. triqs_example:: ising2d_0.cpp
 
-The output is
-
-.. literalinclude:: src/ising2d.output

doc/reference/c++/statistics/ising2d_0.cpp

@@ -0,0 +1,181 @@
+#include <triqs/mc_tools/random_generator.hpp>
+#include <triqs/mc_tools/mc_generic.hpp>
+#include <triqs/utility/callbacks.hpp>
+#include <triqs/arrays.hpp>
+#include <triqs/statistics.hpp>
+#include <vector>
+#include <iostream>
+#include <fstream>
+//#define TRIQS_ARRAYS_ENFORCE_BOUNDCHECK
+// H = -J \sum_<ij> s_i s_j - h \sum_i s_i 
+// theoretical T_c = 2/log(1+sqrt(2)) for J = 1.0
+using namespace triqs::statistics;
+/**************
+ * config 
+ **************/
+
+struct configuration {
+ // N is the linear size of spin matrix, M the total magnetization,
+ // beta the inverse temperature, J the coupling,
+ // field the magnetic field and energy the energy of the configuration
+ int N, M;
+ double beta, J, field, energy;
+ // the chain of spins: true means "up", false means "down"
+ triqs::arrays::array<bool,2> chain;
+ observable<double> M_stack;
+
+ // constructor
+ configuration(int N_, double beta_, double J_, double field_):
+  N(N_), M(-N*N), beta(beta_), J(J_), field(field_), energy(-J*4*N/2+N*field), chain(N,N) , M_stack(){
+   chain()=false;
+  }
+
+};
+
+/**************
+ * move
+ **************/
+
+// A move flipping a random spin
+struct flip {
+ configuration * config;
+ triqs::mc_tools::random_generator &RNG;
+
+ struct site { int i,j ;};//small struct storing indices of a given site
+ site s;
+ double delta_energy;
+
+ // constructor
+ flip(configuration & config_, triqs::mc_tools::random_generator & RNG_) :
+  config(&config_), RNG(RNG_) {}
+
+ // find the neighbours with periodicity
+ std::vector<site> neighbors(site s, int N){
+  std::vector<site> nns(4);
+  int counter=0;
+  for(int i=-1;i<=1;i++){
+   for(int j=-1;j<=1;j++){
+    if ((i==0) != (j==0)) //xor
+     nns[counter++] = site{(s.i+i)%N, (s.j+j)%N};
+   }
+  }
+  return nns;
+ }
+ double attempt() {
+  // pick a random site
+  int index = RNG(config->N*config->N);
+  s = {index%config->N, index/config->N};
+
+  // compute energy difference from field
+  delta_energy = (config->chain(s.i,s.j) ? 2 : -2) * config->field;
+  auto nns = neighbors(s,config->N); //nearest-neighbors
+  double sum_neighbors=0.0;
+  for(auto & x:nns) sum_neighbors += ((config->chain(x.i,x.j))?1:-1);
+  // compute energy difference from J
+  delta_energy += - sum_neighbors * config->J* (config->chain(s.i,s.j)?-2:2);
+
+  // return Metroplis ratio
+  return std::exp(-config->beta * delta_energy);
+ }
+
+ // if move accepted just flip site and update energy and magnetization
+ double accept() {
+  config->M += (config->chain(s.i,s.j) ? -2 : 2);
+  config->chain(s.i,s.j) = !config->chain(s.i,s.j);
+  config->energy += delta_energy;
+  return 1.0;
+ }
+
+ // nothing to do if the move is rejected
+ void reject() {}
+};
+
+
+/**************
+ * measure
+ **************/
+struct compute_m {
+
+ configuration * config;
+ double Z, M;
+
+ compute_m(configuration & config_) : config(&config_), Z(0), M(0) {}
+
+ // accumulate Z and magnetization
+ void accumulate(int sign) {
+
+  Z += sign;
+  M += config->M;
+  //config->M_stack << double(config->M/(config->N*config->N));
+  config->M_stack << config->M;
+ }
+
+ // get final answer M / (Z*N)
+ void collect_results(boost::mpi::communicator const &c) {
+
+  double sum_Z, sum_M;
+  boost::mpi::reduce(c, Z, sum_Z, std::plus<double>(), 0);
+  boost::mpi::reduce(c, M, sum_M, std::plus<double>(), 0);
+
+  if (c.rank() == 0) {
+   std::cout << "@Beta:\t"<<config->beta<<"\tMagnetization:\t" << sum_M / (sum_Z*(config->N*config->N)) << std::endl ;
+   std::cout << "average_and_error(M) = " << average_and_error(config->M_stack) << std::endl;
+   std::cout << "#Beta:\t"<<config->beta<<"\tAutocorr_time:\t" << autocorrelation_time_from_binning(config->M_stack) << std::endl; 
+   std::ofstream outfile("magnetization_series.dat");
+   for(int i=0;i<config->M_stack.size();i++)
+    outfile << config->M_stack[i] <<std::endl;
+   outfile.close();
+  }
+
+ }
+};
+
+int main(int argc, char* argv[]) {
+
+ // initialize mpi
+ boost::mpi::environment env(argc, argv);
+ boost::mpi::communicator world;
+
+ double H=0.0,B=0.5;
+ int N=20;
+ int nc = 100000;
+ if(argc==4){
+  H = atof(argv[1]);//field
+  B = atof(argv[2]);//inverse temp
+  N = atoi(argv[3]);//size along one dimension
+  nc = 1000000 ;
+ }
+ if (world.rank() == 0) 
+  std::cout << "2D Ising with field = " << H << ", beta = " << B << ", N = " << N << std::endl;
+
+ // Prepare the MC parameters
+ int n_cycles = nc;
+ int length_cycle = 100;
+ int n_warmup_cycles = 100000;
+ std::string random_name = "";
+ int random_seed = 374982 + world.rank() * 273894;
+ int verbosity = (world.rank() == 0 ? 2 : 0);
+
+ // Construct a Monte Carlo loop
+ triqs::mc_tools::mc_generic<double> IsingMC(n_cycles, length_cycle, n_warmup_cycles,
+   random_name, random_seed, verbosity);
+
+ // parameters of the model
+ int length = N;
+ double J = 1.0;
+ double field = H;
+ double beta = B;
+
+ // construct configuration
+ configuration config(length, beta, J, field);
+
+ // add moves and measures
+ IsingMC.add_move(flip(config, IsingMC.rng()), "spin flip");
+ IsingMC.add_measure(compute_m(config), "measure magnetization");
+
+ // Run and collect results
+ IsingMC.start(1.0, triqs::utility::clock_callback(-1));
+ IsingMC.collect_results(world);
+
+ return 0;
+}

doc/reference/python/data_analysis/hdf5/tut_ex1.py

@@ -1,11 +0,0 @@
-from pytriqs.archive import *
-import numpy
-R = HDFArchive('myfile.h5', 'w')    # Opens the file myfile.h5, in read/write mode
-R['mu'] = 1.29
-R.create_group('S')
-S= R['S']
-S['a'] = "a string"
-S['b'] = numpy.array([1,2,3])
-del R,S # closing the files (optional : file is closed when the references to R and subgroup are deleted)
-
-

doc/reference/python/data_analysis/hdf5/tut_ex1.rst

@@ -9,7 +9,17 @@ represents the tree structure of the file in a way similar to a dictionary.
 
 Let us start with a very simple example :download:`[file] <./tut_ex1.py>`:
 
-.. literalinclude:: tut_ex1.py
+.. runblock:: python
+
+   from pytriqs.archive import *
+   import numpy
+   R = HDFArchive('myfile.h5', 'w')    # Opens the file myfile.h5, in read/write mode
+   R['mu'] = 1.29
+   R.create_group('S')
+   S= R['S']
+   S['a'] = "a string"
+   S['b'] = numpy.array([1,2,3])
+   del R,S # closing the files (optional : file is closed when the references to R and subgroup are deleted)
 
 Run this and say ::
  

doc/reference/python/data_analysis/hdf5/tut_ex2.py

@@ -1,16 +0,0 @@
-from pytriqs.archive import HDFArchive
-from pytriqs.gf.local import GfImFreq
-
-# Define a Green function 
-G = GfImFreq ( indices = [1], beta = 10, n_points = 1000)
-      
-# Opens the file myfile.h5, in read/write mode
-R = HDFArchive('myfile.h5', 'w')
-# Store the object G under the name 'g1' and mu
-R['g1'] = G
-R['mu'] = 1.29
-del R # closing the files (optional : file is closed when the R reference is deleted)
-
-
-
-

doc/reference/python/data_analysis/hdf5/tut_ex2.rst

@@ -1,4 +1,3 @@
-
 .. _hdf5_tut_ex2:
 
 
@@ -8,13 +7,34 @@ Example 2 : A Green function
 What about more complex objects ?
 The good news is that **hdf-compliant** objects can be stored easily as well.
 
-We can store a Green function in an hdf5 file :download:`[file] <./tut_ex2.py>`:
+We can store a Green function in an hdf5 file:
 
-.. literalinclude:: tut_ex2.py
+.. runblock:: python
 
-Of course, we can retrieve G as easily :download:`[file] <./tut_ex2b.py>`:
+   from pytriqs.archive import HDFArchive
+   from pytriqs.gf.local import GfImFreq
 
-.. literalinclude:: tut_ex2b.py
+   # Define a Green function 
+   G = GfImFreq ( indices = [1], beta = 10, n_points = 1000)
+         
+   # Opens the file myfile.h5, in read/write mode
+   R = HDFArchive('myfile.h5', 'w')
+   # Store the object G under the name 'g1' and mu
+   R['g1'] = G
+   R['mu'] = 1.29
+   del R # closing the files (optional : file is closed when the R reference is deleted)
+
+Of course, we can retrieve G as easily:
+
+.. runblock:: python
+
+   from pytriqs.archive import HDFArchive
+   from pytriqs.gf.local import GfImFreq
+
+   R = HDFArchive('myfile.h5', 'r')  # Opens the file myfile.h5 in readonly mode 
+   G = R['g1'] # Retrieve the object named g1 in the file as G
+
+   # ... ok now I can work with G
 
 The structure of the HDF file is this time ::
 

doc/reference/python/data_analysis/hdf5/tut_ex2b.py

@@ -1,9 +0,0 @@
-from pytriqs.archive import HDFArchive
-from pytriqs.gf.local import GfImFreq
-
-R = HDFArchive('myfile.h5', 'r')  # Opens the file myfile.h5 in readonly mode 
-G = R['g1'] # Retrieve the object named g1 in the file as G
-
-# ... ok now I can work with G
-
-

doc/reference/python/data_analysis/hdf5/tut_ex3.py

@@ -1,12 +0,0 @@
-from pytriqs.gf.local import GfReFreq
-from pytriqs.gf.local.descriptors import SemiCircular
-from pytriqs.archive import HDFArchive
-import numpy
-
-R = HDFArchive('myfile.h5', 'w') 
-for D in range(1,10,2) :
-    g = GfReFreq(indices = [0], window = (-2.00,2.00), name = "D=%s"%D)
-    g <<=  SemiCircular(half_bandwidth = 0.1*D)
-    R[g.name]= g
-
-

doc/reference/python/data_analysis/hdf5/tut_ex3.rst

@@ -11,7 +11,18 @@ Therefore, one can iterate on its elements.
 Let us imagine that you have stored 5 Green functions in an hdf5 file.
 For example, we can create such a file as :download:`[file] <./tut_ex3.py>`:
 
-.. literalinclude:: tut_ex3.py
+.. runblock:: python
+
+   from pytriqs.gf.local import GfReFreq
+   from pytriqs.gf.local.descriptors import SemiCircular
+   from pytriqs.archive import HDFArchive
+   import numpy
+
+   R = HDFArchive('myfile.h5', 'w') 
+   for D in range(1,10,2) :
+       g = GfReFreq(indices = [0], window = (-2.00,2.00), name = "D=%s"%D)
+       g <<=  SemiCircular(half_bandwidth = 0.1*D)
+       R[g.name]= g
 
 
 Imagine that we want to plot those functions :download:`[file] <./tut_ex3b.py>`:

doc/reference/python/lattice/ex_Hilbert.py

@@ -1,28 +0,0 @@
-from pytriqs.lattice.tight_binding import *
-from pytriqs.dos import HilbertTransform
-from pytriqs.gf.local import GfImFreq
-
-# Define a DOS (here on a square lattice)
-BL = BravaisLattice(Units = [(1,0,0) , (0,1,0) ], orbital_positions= {"" :  (0,0,0)} ) 
-t   = -1.00                # First neighbour Hopping
-tp  =  0.0*t               # Second neighbour Hopping
-hop= {  (1,0)  :  [[ t]],       
-        (-1,0) :  [[ t]],     
-        (0,1)  :  [[ t]],
-        (0,-1) :  [[ t]],
-        (1,1)  :  [[ tp]],
-        (-1,-1):  [[ tp]],
-        (1,-1) :  [[ tp]],
-        (-1,1) :  [[ tp]]}
-
-TB = TightBinding (BL, hop)
-d = dos(TB, n_kpts= 500, n_eps = 101, name = 'dos')[0]
-
-#define a Hilbert transform
-H = HilbertTransform(d)
-
-#fill a Green function
-G = GfImFreq(indices = ['up','down'], beta = 20)
-Sigma0 = GfImFreq(indices = ['up','down'], beta = 20); Sigma0.zero()
-G <<= H(Sigma = Sigma0,mu=0.)
-

doc/reference/python/lattice/hilbert.rst

@@ -1,4 +1,3 @@
-
 .. _hilbert_transform:
 
 .. module:: pytriqs.dos.hilbert_transform
@@ -8,8 +7,35 @@ Hilbert Transform
 =======================================
 TRIQS comes with a Hilbert transform. Let us look at an example:
 
-.. literalinclude:: ex_Hilbert.py
+.. runblock:: python
 
+   from pytriqs.lattice.tight_binding import *
+   from pytriqs.dos import HilbertTransform
+   from pytriqs.gf.local import GfImFreq
+
+   # Define a DOS (here on a square lattice)
+   BL = BravaisLattice(units = [(1,0,0) , (0,1,0) ], orbital_positions= [(0,0,0)] ) 
+   t   = -1.00                # First neighbour Hopping
+   tp  =  0.0*t               # Second neighbour Hopping
+   hop= {  (1,0)  :  [[ t]],       
+           (-1,0) :  [[ t]],     
+           (0,1)  :  [[ t]],
+           (0,-1) :  [[ t]],
+           (1,1)  :  [[ tp]],
+           (-1,-1):  [[ tp]],
+           (1,-1) :  [[ tp]],
+           (-1,1) :  [[ tp]]}
+
+   TB = TightBinding (BL, hop)
+   d = dos(TB, n_kpts= 500, n_eps = 101, name = 'dos')[0]
+
+   #define a Hilbert transform
+   H = HilbertTransform(d)
+
+   #fill a Green function
+   G = GfImFreq(indices = ['up','down'], beta = 20)
+   Sigma0 = GfImFreq(indices = ['up','down'], beta = 20); Sigma0.zero()
+   G <<= H(Sigma = Sigma0,mu=0.)
 
 Given a density of states `d` (here for a tight-binding model), the Hilbert transform `H` is defined is defined in the following way::
 

doc/tutorials/python/aim.py

@@ -1,41 +0,0 @@
-from pytriqs.gf.local import *
-from pytriqs.operators import *
-from pytriqs.applications.impurity_solvers.cthyb_matrix import Solver
-
-# Parameters
-D, V, U = 1.0, 0.2, 4.0
-e_f, beta = -U/2.0, 50
-
-# Construct the impurity solver with the inverse temperature
-# and the structure of the Green's functions
-S = Solver(beta = beta, gf_struct = [ ('up',[1]), ('down',[1]) ])
-
-# Initialize the non-interacting Green's function S.G0
-for spin, g0 in S.G0 :
-    g0 <<= inverse( iOmega_n - e_f - V**2 * Wilson(D) ) 
-
-# Run the solver. The result will be in S.G
-S.solve(H_local = U * N('up',1) * N('down',1),   # Local Hamiltonian 
-        quantum_numbers = {                      # Quantum Numbers 
-          'Nup' : N('up',1),                     # Operators commuting with H_Local
-          'Ndown' : N('down',1) },          
-        n_cycles  = 500000,                      # Number of QMC cycles
-        length_cycle = 200,                      # Length of one cycle 
-        n_warmup_cycles = 10000,                 # Warmup cycles
-        n_legendre = 50,                         # Number of Legendre coefficients
-        random_name = 'mt19937',                 # Name of the random number generator
-        use_segment_picture = True,              # Use the segment picture
-        measured_operators = {                   # Operators to be averaged
-          'Nimp' : N('up',1)+N('down',1) }
-        )
-
-# Save the results in an hdf5 file (only on the master node)
-from pytriqs.archive import HDFArchive
-import pytriqs.utility.mpi as mpi
-
-if mpi.is_master_node():
-  Results = HDFArchive("solution.h5",'w')
-  Results["G"] = S.G
-  Results["Gl"] = S.G_legendre
-  Results["Nimp"] = S.measured_operators_results['Nimp']
-

doc/tutorials/python/ctqmc.rst

@@ -35,8 +35,49 @@ we launch the impurity solver calculations by calling the ``Solve`` method.
 Finally, the resulting interacting Green's function as well as average impurity occupancy 
 is stored in the :ref:`HDF5 format<hdf5_base>`.
 
+.. runblock:: python
+
+   from pytriqs.gf.local import *
+   from pytriqs.operators import *
+   from pytriqs.applications.impurity_solvers.cthyb_matrix import Solver
+
+   # Parameters
+   D, V, U = 1.0, 0.2, 4.0
+   e_f, beta = -U/2.0, 50
+
+   # Construct the impurity solver with the inverse temperature
+   # and the structure of the Green's functions
+   S = Solver(beta = beta, gf_struct = [ ('up',[1]), ('down',[1]) ])
+
+   # Initialize the non-interacting Green's function S.G0
+   for spin, g0 in S.G0 :
+       g0 <<= inverse( iOmega_n - e_f - V**2 * Wilson(D) ) 
+
+   # Run the solver. The result will be in S.G
+   S.solve(H_local = U * N('up',1) * N('down',1),   # Local Hamiltonian 
+           quantum_numbers = {                      # Quantum Numbers 
+             'Nup' : N('up',1),                     # Operators commuting with H_Local
+             'Ndown' : N('down',1) },          
+           n_cycles  = 500000,                      # Number of QMC cycles
+           length_cycle = 200,                      # Length of one cycle 
+           n_warmup_cycles = 10000,                 # Warmup cycles
+           n_legendre = 50,                         # Number of Legendre coefficients
+           random_name = 'mt19937',                 # Name of the random number generator
+           use_segment_picture = True,              # Use the segment picture
+           measured_operators = {                   # Operators to be averaged
+             'Nimp' : N('up',1)+N('down',1) }
+           )
+
+   # Save the results in an hdf5 file (only on the master node)
+   from pytriqs.archive import HDFArchive
+   import pytriqs.utility.mpi as mpi
+
+   if mpi.is_master_node():
+     Results = HDFArchive("solution.h5",'w')
+     Results["G"] = S.G
+     Results["Gl"] = S.G_legendre
+     Results["Nimp"] = S.measured_operators_results['Nimp']
 
-.. literalinclude:: ./aim.py
 
 The result can be then read from the ``HDF5`` file and plotted using the ``oplot`` function:
 

doc/tutorials/python/dmft.rst

@@ -20,6 +20,52 @@ the previous single-impurity example to the case of a bath with semi-circular de
 Here is a complete program doing this plain-vanilla DMFT  on a half-filled one-band Bethe lattice:
 
 
-.. literalinclude:: ./dmft.py
+.. runblock:: python
+
+   from pytriqs.gf.local import *
+   from pytriqs.operators import *
+   from pytriqs.archive import *
+   import pytriqs.utility.mpi as mpi
+
+   # Set up a few parameters
+   U = 2.5
+   half_bandwidth = 1.0
+   chemical_potential = U/2.0
+   beta = 100
+   n_loops = 5
+
+   # Construct the CTQMC solver
+   from pytriqs.applications.impurity_solvers.cthyb_matrix import Solver
+   S = Solver(beta = beta, gf_struct = [ ('up',[1]), ('down',[1]) ])
+
+   # Initalize the Green's function to a semi circular
+   S.G <<= SemiCircular(half_bandwidth)
+
+   # Now do the DMFT loop
+   for IterationNumber in range(n_loops):
+
+       # Compute S.G0 with the self-consistency condition while imposing paramagnetism
+       g = 0.5 * ( S.G['up'] + S.G['down'] )
+       for name, g0block in S.G0:
+           g0block <<= inverse( iOmega_n + chemical_potential - (half_bandwidth/2.0)**2  * g )
+
+       # Run the solver
+       S.solve(H_local = U * N('up',1) * N('down',1),                          # Local Hamiltonian
+               quantum_numbers = { 'Nup' : N('up',1), 'Ndown' : N('down',1) }, # Quantum Numbers (operators commuting with H_Local)
+               n_cycles = 5000,                                                # Number of QMC cycles
+               length_cycle = 200,                                             # Length of a cycle
+               n_warmup_cycles = 1000,                                         # How many warmup cycles
+               n_legendre = 30,                                                # Use 30 Legendre coefficients to represent G(tau)
+               random_name = "mt19937",                                        # Use the Mersenne Twister 19937 random generator
+               use_segment_picture = True)                                     # Here we can use the segment picture
+
+       # Some intermediate saves
+       if mpi.is_master_node():
+         R = HDFArchive("single_site_bethe.h5")
+         R["G-%s"%IterationNumber] = S.G
+         del R
+
+       # Here we would usually write some convergence test
+       # if Converged : break
 
 

doc/tutorials/python/ipt.rst

@@ -36,13 +36,47 @@ perturbation theory for a symmetric single-band Anderson model.
 All Green's functions in the calculations have just one index because
 *up* and *down* components are the same.
 
-.. literalinclude:: ipt_solver.py
+.. runblock:: python
+
+   from pytriqs.gf.local import *
+
+   class IPTSolver:
+
+       def __init__(self, **params):
+
+           self.U = params['U']
+           self.beta = params['beta']
+
+           # Matsubara frequency
+           self.g = GfImFreq(indices=[0], beta=self.beta)
+           self.g0 = self.g.copy()
+           self.sigma = self.g.copy()
+
+           # Imaginary time
+           self.g0t = GfImTime(indices=[0], beta = self.beta)
+           self.sigmat = self.g0t.copy()
+
+       def solve(self):
+
+           self.g0t <<= InverseFourier(self.g0)
+           self.sigmat <<= (self.U**2) * self.g0t * self.g0t * self.g0t
+           self.sigma <<= Fourier(self.sigmat)
+
+           # Dyson equation to get G
+           self.g <<= inverse(inverse(self.g0) - self.sigma)
+
 
 Visualization of a Mott transition
 ----------------------------------
 
 We can now use this solver to run DMFT calculations and scan a range of
-values of :math:`U`. At every iteration the resulting data is plotted
+values of :math:`U`. 
+
+.. plot:: tutorials/python/ipt_full.py
+   :include-source:
+   :scale: 70
+
+Alternatively, in this :download:`script <./ipt_dmft.py>`, at every iteration the resulting data is plotted
 and saved into PNG files using the :ref:`TRIQS matplotlib interface<plotting>`.
 Not that :math:`G(i\omega_n)` is analytically continued to the real axis using
 :ref:`Padé approximant<GfReFreq>`.
@@ -50,11 +84,6 @@ Not that :math:`G(i\omega_n)` is analytically continued to the real axis using
 At the end of the script an external utility `convert` is invoked to join the
 DOS plots into a single animated GIF file which illustrates how a metallic
 solution evolves towards an insulator.
-
-.. literalinclude:: ipt_dmft.py
-
-.. only:: html
-
 The result of this script is the following animated gif:
 
   .. image:: mott.gif

doc/tutorials/python/ipt_full.py

@@ -0,0 +1,63 @@
+from pytriqs.gf.local import *
+from pytriqs.plot.mpl_interface import *
+from numpy import *
+import os
+
+class IPTSolver:
+
+    def __init__(self, **params):
+
+        self.U = params['U']
+        self.beta = params['beta']
+
+        # Matsubara frequency
+        self.g = GfImFreq(indices=[0], beta=self.beta)
+        self.g0 = self.g.copy()
+        self.sigma = self.g.copy()
+
+        # Imaginary time
+        self.g0t = GfImTime(indices=[0], beta = self.beta)
+        self.sigmat = self.g0t.copy()
+
+    def solve(self):
+
+        self.g0t <<= InverseFourier(self.g0)
+        self.sigmat <<= (self.U**2) * self.g0t * self.g0t * self.g0t
+        self.sigma <<= Fourier(self.sigmat)
+
+        # Dyson equation to get G
+        self.g <<= inverse(inverse(self.g0) - self.sigma)
+
+# Parameters
+t = 0.5
+beta = 40
+n_loops = 20
+dos_files = []
+
+# Prepare the plot
+plt.figure(figsize=(6,6))
+plt.title("Local DOS, IPT, Bethe lattice, $\\beta=%.2f$"%(beta))
+
+# Main loop over U
+Umax=4.05
+Umin=0.0
+for U in arange(Umin, Umax, 0.51):
+
+    # Construct the IPT solver and set initial G
+    S = IPTSolver(U = U, beta = beta)
+    S.g <<= SemiCircular(2*t)
+
+    # Do the DMFT loop
+    for i in range(n_loops):
+        S.g0 <<= inverse( iOmega_n - t**2 * S.g )
+        S.solve()
+
+    # Get the real-axis with Pade approximants
+    greal = GfReFreq(indices = [1], window = (-4.0,4.0), n_points = 400)
+    greal.set_from_pade(S.g, 201, 0.0)
+
+    r=(U-Umin)/(Umax-Umin) #for color
+    oplot((-1/pi*greal).imag, lw=3,RI='S', color=(r,1.-r,1.-r), label = "U=%1.1f"%U)
+plt.xlim(-4,4)
+plt.ylim(0,0.7)
+plt.ylabel("$A(\omega)$");

foreignlibs/autorun/autorun.py

@@ -74,10 +74,11 @@ class RunBlock(Directive):
         stdout,stderr = proc.communicate(code)
 
         # Process output 
+        out =''
         if stdout:
-            out = ''.join(stdout).decode(output_encoding)
+            out += ''.join(stdout).decode(output_encoding)
         if stderr:
-            out = ''.join(stderr).decode(output_encoding)
+            out += ''.join(stderr).decode(output_encoding)
         
         # Get the original code with prefixes
         if show_source: