diff --git a/doc/reference/python/operators/manipulate.py b/doc/reference/python/operators/manipulate.py new file mode 100644 index 00000000..31b71b93 --- /dev/null +++ b/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')) diff --git a/doc/reference/python/operators/operators.rst b/doc/reference/python/operators/operators.rst new file mode 100644 index 00000000..23f84ba5 --- /dev/null +++ b/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: diff --git a/doc/tutorials/python/ipt/ipt.py b/doc/tutorials/python/ipt/ipt.py new file mode 100644 index 00000000..50b0e3c7 --- /dev/null +++ b/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) diff --git a/doc/tutorials/python/ipt/ipt.rst b/doc/tutorials/python/ipt/ipt.rst new file mode 100644 index 00000000..9a1188a8 --- /dev/null +++ b/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`. 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` and the density of states is plotted into +a PNG file using the :ref:`TRIQS matplotlib interface` +(:math:`G(i\omega_n)` is analytically continued to the real axis by the +:ref:`Padé approximant`). + +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) diff --git a/doc/tutorials/python/ipt/mott.gif b/doc/tutorials/python/ipt/mott.gif new file mode 100644 index 00000000..72a0071c Binary files /dev/null and b/doc/tutorials/python/ipt/mott.gif differ diff --git a/doc/tutorials/python/ipt/mott.py b/doc/tutorials/python/ipt/mott.py new file mode 100644 index 00000000..0a3512cf --- /dev/null +++ b/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)