Add documentation about operators and IPT

new file:   doc/reference/python/operators/
  new file:   doc/tutorials/python/ipt/

doc/reference/python/operators/manipulate.py

@@ -0,0 +1,8 @@
+from pytriqs.operators import *
+
+H = C('up',1) * Cdag('up',2) + C('up',2) * Cdag('up',1)
+print H
+print H - H.dagger()
+
+print anti_commutator(C('up'),Cdag('up'))
+print anti_commutator(C('up'),0.5*Cdag('down'))

doc/reference/python/operators/operators.rst

@@ -0,0 +1,26 @@
+.. _operators:
+
+
+The Operator class
+===================
+
+The TRIQS solvers need to know several operators in order to solve the impurity
+problem. For example, they must know what the local Hamiltonian is, but also its
+quantum numbers (that can be used to improve the speed), possibly some
+operators to be averaged, aso.
+
+In order to deal with these objects, TRIQS provides a class that allows to
+manipulate operators.
+
+A simple example
+-----------------
+
+.. literalinclude:: manipulate.py
+
+
+
+Complete reference
+------------------------
+
+.. autoclass:: pytriqs.operators.Operator
+  :members:

doc/tutorials/python/ipt/ipt.py

@@ -0,0 +1,59 @@
+import numpy
+
+from pytriqs.gf.local import *
+
+class Solver:
+    """A simple IPT solver for the symmetric one band Anderson model"""
+
+    def __init__(self, **params):
+        self.name = 'Iterated Perturbation Theory'
+
+        self.U = params['U']
+        self.beta = params['beta']
+
+        # Only one block in GFs
+        g = GfImFreq(indices =[0], beta =self.beta, name ='0')
+        self.G = BlockGf(name_list=('0',), block_list=(g,))
+        self.G0 = self.G.copy()
+
+    def Solve(self):
+
+        # Imaginary time representation of G_0
+        g0t = GfImTime(indices =[0], beta =self.beta, name ='0')
+        G0t = BlockGf(name_list=('0',), block_list=(g0t,))
+        G0t['0'].set_from_inverse_fourier(self.G0['0'])
+
+        # IPT expression for the self-energy (particle-holy symmetric case is implied)
+        Sigmat = G0t.copy()
+        Sigmat['0'] <<= (self.U**2)*G0t['0']*G0t['0']*G0t['0']
+        
+        self.Sigma = self.G0.copy()
+        self.Sigma['0'].set_from_fourier(Sigmat['0'])
+
+        # New impurity GF from the Dyson's equation
+        self.G <<= self.G0*inverse(1.0 - self.Sigma*self.G0)
+
+S = 0
+
+def run(**params):
+    """IPT loop"""
+
+    # Read input parameters
+    U = params['U']
+    beta = params['beta']
+    N_loops = params['N_loops']
+    Initial_G0 = params['Initial_G0']
+    Self_Consistency = params['Self_Consistency']
+
+    global S
+    # Create a new IPT solver object
+    S = Solver(U=U, beta=beta)
+    # Initialize the bare GF using the function passed in through Initial_G0 parameter
+    Initial_G0(S.G0)
+
+    # DMFT iterations
+    for IterationNumber in range(N_loops):
+        print "DMFT iteration number ", IterationNumber
+
+        S.Solve()
+        Self_Consistency(S.G0,S.G)

doc/tutorials/python/ipt/ipt.rst

@@ -0,0 +1,76 @@
+.. _ipt:
+
+Iterated perturbation theory: an extended DMFT example
+========================================================
+
+Introduction
+------------
+
+The iterated perturbation theory (IPT) was one of the first methods used to solve the
+DMFT equations [#ipt1]_. In spite of its simplistic nature, IPT gives a qualitatively
+correct description of a Mott metal-insulator transition in the Hubbard model on
+infinite-dimensional lattices (on the quantitative level it tends to underestimate
+correlations though). In IPT one iteratively solves the DMFT equations using the
+second-order perturbation theory in Hubbard interaction :math:`U` to approximate
+the impurity self-energy. For the particle-hole symmetric case it reads
+
+.. math::
+
+    \Sigma(i\omega_n) \approx \frac{U}{2} +
+        U^2 \int_0^\beta d\tau e^{i\omega_n\tau} G_0(\tau)^3
+
+A Hartree-Fock contribution :math:`U/2` in the self-energy cancels with a term
+from :math:`G_0(i\omega_n)^{-1}` when the functions are substituted into the
+Dyson's equation. For this reason this contribution is usually omitted from
+both functions.
+
+The success of IPT is caused by the fact that it becomes exact not only in the
+weak coupling limit (by construction), but also reproduces an atomic-limit
+expression for :math:`\Sigma(i\omega_n)` as :math:`U` grows large [#ipt2]_.
+
+Our sample implementation of IPT includes two files: `ipt.py` and
+`mott.py`.
+
+IPT solver and self-consistency loop
+------------------------------------
+
+The file `ipt.py` implements the weak coupling perturbation theory for a 
+symmetric single-band Anderson model (`Solver` class) and obeys
+:ref:`the solver concept<solver_concept>`. It also runs a DMFT loop with this
+solver and with a self-consistency condition provided from outside (function `run`).
+
+All Green's functions in the calculations consist of one Matsubara block
+(there is no need for spin indices, since only paramagnetic solutions are sought).
+
+.. literalinclude:: ipt.py
+
+Visualization of a Mott transition
+----------------------------------
+
+In `mott.py` the IPT module is used to run DMFT many times and scan a range of 
+values of :math:`U`. On every run the resulting data are saved in 
+an :ref:`HDF5 archive<hdf5_base>` and the density of states is plotted into
+a PNG file using the :ref:`TRIQS matplotlib interface<plotting>`
+(:math:`G(i\omega_n)` is analytically continued to the real axis by the
+:ref:`Padé approximant<GfReFreq>`). 
+
+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:: mott.py
+
+The result of this script is an animated gif:
+
+.. image:: mott.gif
+   :width: 700
+   :align: center
+
+
+Journal references
+------------------
+
+.. [#ipt1] A. Georges and G. Kotliar,
+           Phys. Rev. B 45, 6479–6483 (1992).
+.. [#ipt2] X. Y. Zhang, M. J. Rozenberg, and G. Kotliar,
+           Phys. Rev. Lett. 70, 1666–1669 (1993)

doc/tutorials/python/ipt/mott.py

@@ -0,0 +1,81 @@
+# Visualization of the Mott transition
+
+from math import *
+import os
+import numpy
+
+from pytriqs.gf.local import *
+from pytriqs.gf.local import Omega, SemiCircular, inverse
+from pytriqs.archive import *
+from pytriqs.plot.mpl_interface import oplot
+import matplotlib.pyplot as plt
+
+import ipt
+
+beta = 40      # Inverse temperature
+U = numpy.arange(0,4.05,0.1)    # Range of interaction constants to scan
+t = 0.5     # Scaled hopping constant on the Bethe lattice
+
+N_loops = 20    # Number of DMFT loops
+
+# Pade-related
+DOSMesh = numpy.arange(-4,4,0.02)   # Mesh to plot densities of states
+eta = 0.00     # Imaginary frequency offset to use with Pade
+Pade_L = 201    # Number of Matsubara frequencies to use with Pade
+
+# Bare Green's function of the Bethe lattice (semicircle) 
+def Initial_G0(G0):
+    G0 <<= SemiCircular(2*t)
+
+# Self-consistency condition
+def Self_Consistency(G0,G):
+    G0['0'] <<= inverse(Omega - (t**2)*G['0'])
+
+# Save results to an HDF5-archive
+ar = HDFArchive('Mott.h5','w')
+ar['beta'] = beta
+ar['t'] = t
+    
+# List of files with DOS figures
+DOS_files=[]
+    
+# Scan over values of U
+for u in U:
+    print "Running DMFT calculation for beta=%.2f, U=%.2f, t=%.2f" % (beta,u,t)
+    ipt.run(N_loops = N_loops, beta=beta, U=u, Initial_G0=Initial_G0, Self_Consistency=Self_Consistency)
+    
+    # The resulting local GF on the real axis
+    g_real = GfReFreq(indices = [0], beta = beta, mesh_array = DOSMesh, name = '0')
+    G_real = BlockGf(name_list = ('0',), block_list = (g_real,), make_copies = True)
+    
+    # Analytic continuation with Pade
+    G_real['0'].set_from_pade(ipt.S.G['0'], N_Matsubara_Frequencies=Pade_L, Freq_Offset=eta)
+
+    # Save data to the archive
+    ar['U' + str(u)] = {'G0': ipt.S.G0, 'G': ipt.S.G, 'Sigma': ipt.S.Sigma, 'G_real':G_real}
+   
+    # Plot the DOS
+    fig = plt.figure()
+    oplot(G_real['0'][0,0], RI='S', name="DOS", figure = fig)
+    
+    # Adjust 'y' axis limits accordingly to the Luttinger sum rule 
+    fig.axes[0].set_ylim(0,1/pi/t*1.1)
+    
+    # Set title of the plot
+    fig_title = "Local DOS, IPT, Bethe lattice, $\\beta=%.2f$, $U=%.2f$" % (beta,u)
+    plt.title(fig_title)
+
+    # Save the figure as a PNG file
+    DOS_file = "DOS_beta%.2fU%.2f.png" % (beta,u)
+    fig.savefig(DOS_file, format="png", transparent=False)
+    DOS_files.append(DOS_file)
+    plt.close(fig)
+
+# Create an animated GIF
+# (you need to have 'convert' utility installed; it is a part of ImageMagick suite) 
+convert_cmd = "convert -delay 25 -loop 0"
+convert_cmd += " " + ' '.join(DOS_files)
+convert_cmd += " " + "DOS.gif"
+
+print "Creating an animated DOS plot..."
+os.system(convert_cmd)