HUGE UPDATE: quite deep update/maquillage of the code following the creation of a new mixgf_dmft and upm_mixgf_dmft branches, where I introduced new features and practices. The code is successfully tested against U=4,E=1,V=1 case.
Restarting stategy: read Sigma_less/gtr file AND nk file.
All other guess (G0local, G0local+electric field) have been commented out
(not erased to be restored if needed)
Important Update: collision integrals are evaluated using different and possibly more precise
methods: rectangular, trapezoid and simpson rule. A flag int_method (default trapz) is introduced.
Result are as follow. trapz and simps gives very similar results, but the two are somewhat
different from the result with rect. rule (expected given the increased precision).
No update: just updated inputFILE to do some calculations for Stark-Ladder
UPDATED code: small revision, updated with improvements from upm_mixgf_dmft branch. Tested against previous results successfull. Starting point for branching to nomix_devel_hf
Update: using update_WF method #3: direct inversion of the Keldysh GF.
UPDATE: added more bath type: pseudo-gap and gapped bath. Tests running.
update: last changes to evaluate the optical conductivity. not finished.
Some bug fixing due to scifor library update.
rapid commit to debug devel_neq_code_rk4.
updated code, simplified structure.
Commit to branch to new Kadanoff-Baym solver
based on collocation method.
small update to this old version of the code to compile and have a working executable.
Ready to merge with nomix_devel_newcode
update: some changes here. I discovered a bug introduced in the
discretization of the funcional array. working on it.

New version of the code.
Kadanoff-Baym equations for the evaluation of the
local Green's functions are sovled using Runge-Kutta
/Collocation method for Volterra Integro-Differential Equations (VIDE).

WARNING: the code is not yet test. Only the U=0 has been studied.
More tests are required, specially to test self-consistency condition.
update: tested U=2 self-consistency, 8 DMFT-loop with no problems.
checking correctness of the results.

The actual solution is obtained with 2nd order method.
Higher order will be implemented soon (still debugging the test driver).

Notes:
- While updating the code some changes have been implemented such as rename of
the modules to include NEQ_ statement so to avoid overlap with other codes.

- The entire code has been ported to the 1:nstep rather than 0:nstep formalism.

- Equilibrium solver is about to be entirely removed. Module is in local repo.

- The stationary or equilibrium GF have been commented out. I should write
the correct extension to the 1-nstep:nstep-1 basis. Working on it.

TODO:
1. To save roughly half of the memory we can take advantage of the GF symmetries and store the
GF into lower/upper triangular part. This can be copied from the RK4 code in
the *devel_rk4kb_new_code* branch, using RK implicit method.

2. The hamiltonian must be passed as abstract procedure so to make the
KB routine generic w/ respect to Ham

3. Code higher order in weak-coupling perturbation theory solvers.
updated code: the new version had a bug in the Kadanoff-Baym equatuion solution
introduced with the discretization of the functions arguments in the
vide_rkX routines.
I restored the older version with function as arguments and checked that
the code gives exactly the same results in the non-interacting limit as
the test drivers testGFviderk2 and neqGZ codes.

update: testing the U!=0 case.
G.Mazza suggested that the <-component should be calculated only in the
upper triangular regime, which decreases the calculation time.

The possible large error along the diagonal, related to the sum of diagonal in k,
should come from the double errors: one on the initial conditions and one over the dynamics.
The former are in fact obtained from the dynamics along the first line (t`=0).
Can we do better?
just a commit, I am w/ Gabi
debug: keep testing this new version of the code.
I fixed some issues with the evaluation of Sigma function (IPT) from
G0. The gtr component was wrongly evaluated from less and ret component
which generated a bad Sigma function, compromising the entire calculation.

Other fixing in the code in BATH, KADANOFF-BAYM, IPT.

Rational @260913:
I am testing the iteration with and without the interactino and bath self-energy.
Somehow (and here is the bug I think), in presence of BOTH the self-energies
the system does not conserve the number of particles. In passing I should notice that
this bug remind very much to the gapped/pseud-gapped bath issue of non-conserving density.
If U=0/V!=0 OR V=0/U!=0 the density is conserved. Note that in the first case we are
including the effect of coupling quench with the thermostat.
update: removed the Qkret terms proportional to delta function.
After discussino with G.Mazza, we realized that the delta terms should
not be there becuase the solution is obtained for t>t', that is
we are solving in the left open interval.

bug: major (dummy!) bug found and fixed.
The Kernel function Skernel in the KB equations evolution
should include the tstride so be correct.
This is because the Ret. component has the integral going
from t'---> t == ytime:Nstep, while the functions passed to *vide_rk2* routine
as a reference has indices going from 1:nt.
Thus when calling the function you used a kernel shifted back, which
in presence of broken translation invariance (non-equilibrium), i.e.
a non-eq Sigma function different from the Bath, produced wrong results.
Now fixed.

One small little bug remains I guess, because the code misbehave at larger U.
NEW VERSION OF THE CODE.
This version is based on causal evoution of the local function.
All the equations are reduced to integral or integro-differential equations.

A large class/module for the contour GF is written, including mixed gf.

The code is now very short and all included in almost one file in the driver dir.

In this version it keeps using IPT solver. As a starting guess we use G0 (the Weiss Field)
but this requires the user to write a specific junction between input and the solver, basically
reshaping the dmft-loop with respect to my usual implmentation (SC--->IS), now is IS--->SC.

Checking the stability of the code it seems like that the code performs better at larger U if the
matsubara/imaginary time integral are performed with a larger basis. Namely because the self-energy
is very sharp at the boundaries (0,beta). I will implement a version of the code using PQ-basis.

Another driver is in the middle of development to start from a given self-energy Matsubara.
This requires to Fourier transforms to imaginary time, which is only possible if nasty tail C_0 + C_1/iw
are removed and transformed alone.
Checking on this.
CODE WORKING
The new version of the code using causal propagation at 2nd order and
VIE/VIDE solution has been thoroughly tested against the available data
for the quench.
[This version of the code is based on the algorithm provided by
M.Eckstein, P.Werner, N.Tsuji and Co., plus some minor improvements]

The main driver: +neqdmft_bethe_dos was shown to give identical results
as its related code neqdmft_bethe performing algebraic self-consistency on
the bethe lattice. The results are identical to published data.

Tests for the quench in the 2D square lattice are ok.
Test for the electric field effects on the 2D square lattice are going.
UPDATE: the solution at large U are affected by breakdown of the
density conservation! (both IC). This is certainly related to the
poor quality of the IPT solver.
Running tests to investigate the reliability of the code.

updated version: working with correlated initial conditions.
Can not yet solve large U limit, that indeed seems to work starting from
uncorrelated initial conditions.

The solutions (un/correlated IC) reach the same steady state, checked from U=4,8 in 2D
square lattice.
NEW VERSION OF THE CODE.
This version uses causal algorithm + VIE/VIDE solution

The routines are adaption of the previous version of the code
to the method outlined by Eckstein, Tsuji, Werner & Co. and
presented in their sample code, attached to the paper:
1310.5329 // available at http://arxiv.org/abs/1310.5329v1

Merge branch 'devel_dmft'

Conflicts:
	CONTOUR_GF.f90
	FUNX_NEQ.f90
	IPT_NEQ.f90
	KADANOFBAYM.f90
	Makefile
	VARS_GLOBAL.f90
	neqDMFT.f90
	update_G0_nonequilibrium.f90