#if defined HAVE_CONFIG_H
#include "config.h"
#endif

!{\src2tex{textfont=tt}}
!!****m* ABINIT/m_Ctqmc
!! NAME
!!  m_Ctqmc
!! 
!! FUNCTION 
!!  Manage and drive all the CTQMC
!!  Should not be used if you don't know what you do
!!  Please use CtqmcInterface
!!
!! COPYRIGHT
!!  Copyright (C) 2013-2014 ABINIT group (J. Bieder)
!!  This file is distributed under the terms of the
!!  GNU General Public License, see ~abinit/COPYING
!!  or http://www.gnu.org/copyleft/gpl.txt .
!!
!! NOTES
!!
!! PARENTS
!!  Will be filled automatically by the parent script
!!
!! CHILDREN
!!  Will be filled automatically by the parent script
!!
!! SOURCE

#include "defs.h"

MODULE m_Ctqmc

USE m_Global
USE m_GreenHyb
USE m_BathOperator
USE m_ImpurityOperator
USE m_Stat
USE m_FFTHyb
USE m_OurRng
USE m_Vector
#ifdef HAVE_MPI2
USE mpi
#endif

IMPLICIT NONE

!!***

PRIVATE

INTEGER, PARAMETER :: CTQMC_SLICE1 = 100
! Coupe Sweeps en 100
INTEGER, PARAMETER :: CTQMC_SLICE2 = 100
! Coupe modNoise1 en 100
INTEGER, PARAMETER :: CTQMC_SEGME =  1
INTEGER, PARAMETER :: CTQMC_ANTIS = -2
INTEGER, PARAMETER :: CTQMC_ADDED =  3  
INTEGER, PARAMETER :: CTQMC_REMOV =  4
INTEGER, PARAMETER :: CTQMC_DETSI =  5


!!****t* m_Ctqmc/Ctqmc
!! NAME
!!  Ctqmc
!!
!! FUNCTION
!!  This structured datatype contains the necessary data
!!
!! COPYRIGHT
!!  Copyright (C) 2013-2014 ABINIT group (J. Bieder)
!!  This file is distributed under the terms of the
!!  GNU General Public License, see ~abinit/COPYING
!!  or http://www.gnu.org/copyleft/gpl.txt .
!!
!! SOURCE

TYPE, PUBLIC :: Ctqmc

  LOGICAL _PRIVATE :: init = .FALSE.
! Flag: is MC initialized

  LOGICAL _PRIVATE :: set  = .FALSE.
! Flag: ??

  LOGICAL _PRIVATE :: setU = .FALSE.
! Flag: is U Set ?

  LOGICAL _PRIVATE :: inF  = .FALSE.
! Flag: is hybridization fct in input ?

  LOGICAL _PRIVATE :: done = .FALSE.
! Flag: is MC terminated ?

  LOGICAL _PRIVATE :: para = .FALSE.
! Flag:  do we have parameters in input

  LOGICAL _PRIVATE :: have_MPI = .FALSE.
! Flag: 

  INTEGER _PRIVATE :: opt_movie = 0
!

  INTEGER _PRIVATE :: opt_analysis = 0
! correlations 

  INTEGER _PRIVATE :: opt_check = 0
! various check 0
! various check 1 impurity
! various check 2 bath
! various check 3 both

  INTEGER _PRIVATE :: opt_order = 0
! nb of segments max for analysis

  INTEGER _PRIVATE :: opt_noise = 0
! compute noise

  INTEGER _PRIVATE :: opt_spectra = 0
! markov chain FT (correlation time)

  INTEGER _PRIVATE :: opt_levels = 0
! do we have energy levels

  INTEGER _PRIVATE :: flavors
!

  INTEGER _PRIVATE :: measurements
! nb of measure in the MC

  INTEGER _PRIVATE :: samples
! nb of L points

  INTEGER(8) _PRIVATE :: seed
!

  INTEGER _PRIVATE :: sweeps
!

  INTEGER _PRIVATE :: thermalization
!

  INTEGER _PRIVATE :: ostream
! output file

  INTEGER _PRIVATE :: istream
! input file

  INTEGER _PRIVATE :: modNoise1
! measure the noise each modNoise1

  INTEGER _PRIVATE :: modNoise2
! measure the noise each modNoise2

  INTEGER _PRIVATE :: activeFlavor
! orbital on which one do sth now

  INTEGER, DIMENSION(1:2) _PRIVATE :: modGlobalMove
! 1: gloabl move each modglobalmove(1)
! 2: we have done modglobalmove(2) for two different orbitals.

  INTEGER _PRIVATE :: Wmax
! Max freq for FT

  DOUBLE PRECISION, DIMENSION(1:6) _PRIVATE :: stats
! to now how many negative determinant, antisegments,seeme.e.twfs...j

  DOUBLE PRECISION _PRIVATE :: swap
! nb of successfull GM

  INTEGER _PRIVATE :: MY_COMM
! 

  INTEGER _PRIVATE :: rank
!

  INTEGER _PRIVATE :: size
! size of MY_COMM

  DOUBLE PRECISION _PRIVATE :: runTime ! time for the run routine
!  

  DOUBLE PRECISION _PRIVATE :: beta
!

  DOUBLE PRECISION _PRIVATE :: U

  DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:) _PRIVATE :: mu
! levels

  TYPE(GreenHyb)  , ALLOCATABLE, DIMENSION(:    ) _PRIVATE :: Greens 

  DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:  ) _PRIVATE :: measN 
! measure of occupations (3or4,flavor) 

  DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:  ) _PRIVATE :: measDE
!  (flavor,flavor) double occupancies
!  (1,1): total energy of correlation.

  DOUBLE PRECISION _PRIVATE :: a_Noise
! Noise a exp (-bx) for the  noise

  DOUBLE PRECISION _PRIVATE :: b_Noise
! Noise a exp (-bx) for the  noise

  DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:,:) _PRIVATE :: abNoiseG   !(ab,tau,flavor)
! Noise but for G

  TYPE(Vector)             , DIMENSION(1:2) _PRIVATE :: measNoise 
  TYPE(Vector), ALLOCATABLE, DIMENSION(:,:,:) _PRIVATE :: measNoiseG       !(tau,flavor,mod) 
! accumulate each value relataed to measurenoise 1 2

  DOUBLE PRECISION _PRIVATE                            :: inv_dt
! 1/(beta/L)

  DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:  ) _PRIVATE :: measPerturbation 
! opt_order,nflavor

  DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:,:) _PRIVATE :: measCorrelation 
! segment,antisegment,nflavor,nflavor

  DOUBLE PRECISION _PRIVATE :: errorImpurity
! check 

  DOUBLE PRECISION _PRIVATE :: errorBath
! for check

  TYPE(BathOperator) _PRIVATE              :: Bath

  TYPE(ImpurityOperator) _PRIVATE          :: Impurity

  DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:) _PRIVATE :: density

END TYPE Ctqmc
!!***

PUBLIC  :: Ctqmc_init
PRIVATE :: Ctqmc_setParameters
PUBLIC  :: Ctqmc_setSweeps
PUBLIC  :: Ctqmc_setSeed
PRIVATE :: Ctqmc_allocateAll
PRIVATE :: Ctqmc_allocateOpt
PUBLIC  :: Ctqmc_setG0wFile
PUBLIC  :: Ctqmc_setG0wTab
PUBLIC  :: Ctqmc_setU
PRIVATE :: Ctqmc_clear
PUBLIC  :: Ctqmc_reset
PUBLIC  :: Ctqmc_setMu
PRIVATE :: Ctqmc_computeF
PUBLIC  :: Ctqmc_run
PRIVATE :: Ctqmc_loop
PRIVATE :: Ctqmc_tryAddRemove
PRIVATE :: Ctqmc_trySwap
PRIVATE :: Ctqmc_measN
PRIVATE :: Ctqmc_measCorrelation
PRIVATE :: Ctqmc_measPerturbation
PUBLIC  :: Ctqmc_getResult
PUBLIC  :: Ctqmc_symmetrizeGreen
PUBLIC  :: Ctqmc_getGreen
PUBLIC  :: Ctqmc_getD
PUBLIC  :: Ctqmc_getE
PUBLIC  :: Ctqmc_printAll
PUBLIC  :: Ctqmc_printQMC
PUBLIC  :: Ctqmc_printGreen
PUBLIC  :: Ctqmc_printD
PUBLIC  :: Ctqmc_printE
PUBLIC  :: Ctqmc_printPerturbation
PUBLIC  :: Ctqmc_printCorrelation
PUBLIC  :: Ctqmc_printSpectra
PUBLIC  :: Ctqmc_destroy

CONTAINS
!!***

!!****f* ABINIT/m_Ctqmc/Ctqmc_init
!! NAME
!!  Ctqmc_init
!!
!! FUNCTION
!!  Initialize the type Ctqmc
!!  Allocate all the non optional variables
!!
!! COPYRIGHT
!!  Copyright (C) 2013-2014 ABINIT group (J. Bieder)
!!  This file is distributed under the terms of the
!!  GNU General Public License, see ~abinit/COPYING
!!  or http://www.gnu.org/copyleft/gpl.txt .
!!
!! INPUTS
!!  op=ctqmc
!!  ostream=where to write
!!  istream=where to read the input parameters if so
!!  bFile=do we read in istream ?
!!  MY_COMM=mpi communicator for the CTQMC
!!  iBuffer=input parameters if bFile is false
!!
!! OUTPUT
!!
!! SIDE EFFECTS
!!
!! NOTES
!!
!! PARENTS
!!  Will be filled automatically by the parent script
!!
!! CHILDREN
!!  Will be filled automatically by the parent script
!!
!! SOURCE

SUBROUTINE Ctqmc_init(op, ostream, istream, bFile, MY_COMM, iBuffer)


!This section has been created automatically by the script Abilint (TD).
!Do not modify the following lines by hand.
#undef ABI_FUNC
#define ABI_FUNC 'Ctqmc_init'
!End of the abilint section


#ifdef HAVE_MPI1
include 'mpif.h'
#endif
!Arguments ------------------------------------
  TYPE(Ctqmc), INTENT(INOUT)                      :: op
  INTEGER  , INTENT(IN   )                      :: ostream
  INTEGER  , INTENT(IN   )                      :: istream
  LOGICAL  , INTENT(IN   )                      :: bFile
  DOUBLE PRECISION, DIMENSION(1:9), OPTIONAL, INTENT(IN) :: iBuffer
  INTEGER  , OPTIONAL, INTENT(IN   )                      :: MY_COMM
!Local variables ------------------------------
#ifdef HAVE_MPI
  INTEGER                                       :: ierr
#endif
  INTEGER                                       :: iflavor
#ifdef __GFORTRAN__
!  INTEGER                                       :: pid
!  CHARACTER(LEN=5)                              :: Cpid
!
#endif
  DOUBLE PRECISION, DIMENSION(1:9)             :: buffer

  op%ostream = ostream
  op%istream = istream
  
! --- RENICE ---
!#ifdef __GFORTRAN__
!  pid = GetPid()
!  WRITE(Cpid,'(I5)') pid
!  CALL SYSTEM('renice +19 '//TRIM(ADJUSTL(Cpid))//' > /dev/null')
!#endif
!! --- RENICE ---

  IF ( PRESENT(MY_COMM)) THEN
#ifdef HAVE_MPI
    op%have_MPI = .TRUE.
    op%MY_COMM = MY_COMM
    CALL MPI_Comm_rank(op%MY_COMM, op%rank, ierr)
    CALL MPI_Comm_size(op%MY_COMM, op%size, ierr)
#else
    CALL WARN("MPI is not used                                    ")
    op%have_MPI = .FALSE.
    op%MY_COMM = -1
    op%rank = 0
    op%size = 1
#endif
  ELSE
    op%have_MPI = .FALSE.
    op%MY_COMM = -1
    op%rank = 0
    op%size = 1
  END IF

  !IF ( op%rank .EQ. 0 ) THEN
!  WRITE(ostream,'(A20)') 'Job reniced with +19'
    !CALL FLUSH(ostream)
  !END IF

  IF ( bFile .EQV. .TRUE. ) THEN
    IF ( op%rank .EQ. 0 ) THEN

      READ(istream,*) buffer(1) !iseed
      READ(istream,*) buffer(2) !op%sweeps
      READ(istream,*) buffer(3) !op%thermalization
      READ(istream,*) buffer(4) !op%measurements
      READ(istream,*) buffer(5) !op%flavors
      READ(istream,*) buffer(6) !op%samples
      READ(istream,*) buffer(7) !op%beta
      READ(istream,*) buffer(8) !U
      READ(istream,*) buffer(9) !iTech
      !READ(istream,*) buffer(9) !Wmax
!#ifdef CTCtqmc_ANALYSIS
      !READ(istream,*) buffer(10) !order
!#endif
    END IF

#ifdef HAVE_MPI
    IF ( op%have_MPI .EQV. .TRUE. ) &
      CALL MPI_Bcast(buffer, 9, MPI_DOUBLE_PRECISION, 0,    &
                   op%MY_COMM, ierr)
#endif
  ELSE IF ( PRESENT(iBuffer) ) THEN
    buffer(1:9) = iBuffer(1:9)
  ELSE
    CALL ERROR("Ctqmc_init : No input parameters                    ")
  END IF

  CALL Ctqmc_setParameters(op, buffer)

  CALL Ctqmc_allocateAll(op)

  DO iflavor = 1, op%flavors
      CALL GreenHyb_init(op%Greens(iflavor),op%samples, op%beta, iTech=INT(buffer(9)),MY_COMM=op%MY_COMM)
  END DO


!  op%seg_added    = 0.d0
!  op%anti_added   = 0.d0
!  op%seg_removed  = 0.d0
!  op%anti_removed = 0.d0
!  op%seg_sign     = 0.d0
!  op%anti_sign    = 0.d0
  op%stats(:)     = 0.d0
  op%swap         = 0.d0
  op%runTime      = 0.d0

  CALL Vector_init(op%measNoise(1),op%sweeps/op%modNoise1)
  CALL Vector_init(op%measNoise(2),(op%sweeps/op%modNoise1+1)*CTQMC_SLICE2)
  !CALL Vector_init(op%measNoise(3),101)
  !CALL Vector_init(op%measNoise(4),101)

  op%set  = op%para .AND. op%inF
  op%done = .FALSE.
  op%init = .TRUE.

!#ifdef CTCtqmc_CHECK
  op%errorImpurity = 0.d0
  op%errorBath     = 0.d0
!#endif
END SUBROUTINE Ctqmc_init
!!***

!!****f* ABINIT/m_Ctqmc/Ctqmc_setParameters
!! NAME
!!  Ctqmc_setParameters
!!
!! FUNCTION
!!  set all parameters and operators
!!
!! COPYRIGHT
!!  Copyright (C) 2013-2014 ABINIT group (J. Bieder)
!!  This file is distributed under the terms of the
!!  GNU General Public License, see ~abinit/COPYING
!!  or http://www.gnu.org/copyleft/gpl.txt .
!!
!! INPUTS
!!  op=ctqmc
!!  buffer=input parameters
!!
!! OUTPUT
!!
!! SIDE EFFECTS
!!
!! NOTES
!!
!! PARENTS
!!  Will be filled automatically by the parent script
!!
!! CHILDREN
!!  Will be filled automatically by the parent script
!!
!! SOURCE

SUBROUTINE Ctqmc_setParameters(op,buffer)

!Arguments ------------------------------------

!This section has been created automatically by the script Abilint (TD).
!Do not modify the following lines by hand.
#undef ABI_FUNC
#define ABI_FUNC 'Ctqmc_setParameters'
!End of the abilint section

  TYPE(Ctqmc), INTENT(INOUT)                         :: op
  DOUBLE PRECISION, DIMENSION(1:9), INTENT(IN   ) :: buffer


  op%thermalization = INT(buffer(3)) !op%thermalization
  CALL Ctqmc_setSeed(op,INT(buffer(1)))
  CALL Ctqmc_setSweeps(op,buffer(2))

  op%measurements   = INT(buffer(4)) !op%measurements
  op%flavors        = INT(buffer(5)) !op%flavors
  op%samples        = INT(buffer(6)) !op%samples
  op%beta           = buffer(7)      !op%beta
  op%U              = buffer(8)      !U
!  op%mu             = buffer(9)      !op%mu
  !op%Wmax           = INT(buffer(9)) !Freq
!#ifdef CTCtqmc_ANALYSIS
!  op%order          = INT(buffer(10)) ! order
  op%inv_dt         = op%samples / op%beta
!#endif

  !CALL ImpurityOperator_init(op%Impurity,op%flavors,op%beta, op%samples)
  CALL ImpurityOperator_init(op%Impurity,op%flavors,op%beta)
  IF ( op%U .GE. 0.d0 ) THEN
    CALL ImpurityOperator_computeU(op%Impurity,op%U,0.d0)
    op%setU = .TRUE.
  END IF
!  op%mu = op%mu + op%Impurity%shift_mu

  CALL BathOperator_init(op%Bath, op%flavors, op%samples, op%beta, INT(buffer(9)))

  op%para = .TRUE.

END SUBROUTINE Ctqmc_setParameters
!!***

!!****f* ABINIT/m_Ctqmc/Ctqmc_setSweeps
!! NAME
!!  Ctqmc_setSweeps
!!
!! FUNCTION
!!  set the number of sweeps
!!
!! COPYRIGHT
!!  Copyright (C) 2013-2014 ABINIT group (J. Bieder)
!!  This file is distributed under the terms of the
!!  GNU General Public License, see ~abinit/COPYING
!!  or http://www.gnu.org/copyleft/gpl.txt .
!!
!! INPUTS
!!  op=ctqmc
!!  sweeps=asked sweeps
!!
!! OUTPUT
!!
!! SIDE EFFECTS
!!
!! NOTES
!!
!! PARENTS
!!  Will be filled automatically by the parent script
!!
!! CHILDREN
!!  Will be filled automatically by the parent script
!!
!! SOURCE

SUBROUTINE Ctqmc_setSweeps(op,sweeps)

!Arguments ------------------------------------

!This section has been created automatically by the script Abilint (TD).
!Do not modify the following lines by hand.
#undef ABI_FUNC
#define ABI_FUNC 'Ctqmc_setSweeps'
!End of the abilint section

  TYPE(Ctqmc)         , INTENT(INOUT) :: op
  DOUBLE PRECISION  , INTENT(IN   ) :: sweeps

  op%sweeps = NINT(sweeps / DBLE(op%size))
!  write(6,*) op%sweeps,NINT(sweeps / DBLE(op%size)),ANINT(sweeps/DBLE(op%size))
  IF ( DBLE(op%sweeps) .NE. ANINT(sweeps/DBLE(op%size)) ) &
    CALL ERROR("Ctqmc_setSweeps : sweeps is negative or too big     ")
  IF ( op%sweeps .LT. 2*CTQMC_SLICE1 ) THEN  !202
    CALL WARNALL("Ctqmc_setSweeps : # sweeps automtically changed     ")
    op%sweeps = 2*CTQMC_SLICE1
!  ELSE IF ( op%sweeps .LT. op%thermalization ) THEN
!    CALL WARNALL("Ctqmc_setSweeps : Thermalization > sweeps / cpu -> auto fix")
!    op%sweeps = op%thermalization
  END IF
  IF ( DBLE(NINT(DBLE(op%sweeps)*DBLE(op%size)/DBLE(CTQMC_SLICE1))) .NE.  &
  ANINT(DBLE(op%sweeps)*DBLE(op%size)/DBLE(CTQMC_SLICE1)) ) THEN
    op%modNoise1 = op%sweeps
  ELSE
    op%modNoise1    = MIN(op%sweeps,INT(DBLE(op%sweeps)*DBLE(op%size) / DBLE(CTQMC_SLICE1))) !101
  END IF
  op%modNoise2    = MAX(op%modNoise1 / CTQMC_SLICE2, 1)   ! 100
!  op%modGlobalMove(1) = op%thermalization / 10 + 1
!  op%modGlobalMove(2) = 0

END SUBROUTINE Ctqmc_setSweeps
!!***

!!****f* ABINIT/m_Ctqmc/Ctqmc_setSeed
!! NAME
!!  Ctqmc_setSeed
!!
!! FUNCTION
!!  initialize random number generator
!!
!! COPYRIGHT
!!  Copyright (C) 2013-2014 ABINIT group (J. Bieder)
!!  This file is distributed under the terms of the
!!  GNU General Public License, see ~abinit/COPYING
!!  or http://www.gnu.org/copyleft/gpl.txt .
!!
!! INPUTS
!!  op=ctqmc
!!  iseed=seed from imput
!!
!! OUTPUT
!!
!! SIDE EFFECTS
!!
!! NOTES
!!
!! PARENTS
!!  Will be filled automatically by the parent script
!!
!! CHILDREN
!!  Will be filled automatically by the parent script
!!
!! SOURCE

SUBROUTINE Ctqmc_setSeed(op,iseed)

!Arguments ------------------------------------

!This section has been created automatically by the script Abilint (TD).
!Do not modify the following lines by hand.
#undef ABI_FUNC
#define ABI_FUNC 'Ctqmc_setSeed'
!End of the abilint section

  TYPE(Ctqmc), INTENT(INOUT)           :: op
  INTEGER  , INTENT(IN   )           :: iseed
!Local variables ------------------------------
  !INTEGER                            :: n
  !INTEGER                            :: i
  !INTEGER, DIMENSION(:), ALLOCATABLE :: seed


  !CALL RANDOM_SEED(size = n)
  !MALLOC(seed,(n))
  !seed =  iseed + (/ (i - 1, i = 1, n) /)

  !CALL RANDOM_SEED(PUT = seed+op%rank)

  !FREE(seed)

  op%seed=INT(iseed+op%rank,8)

END SUBROUTINE Ctqmc_setSeed
!!***

!!****f* ABINIT/m_Ctqmc/Ctqmc_allocateAll
!! NAME
!!  Ctqmc_allocateAll
!!
!! FUNCTION
!!  Allocate all non option varibales
!!
!! COPYRIGHT
!!  Copyright (C) 2013-2014 ABINIT group (J. Bieder)
!!  This file is distributed under the terms of the
!!  GNU General Public License, see ~abinit/COPYING
!!  or http://www.gnu.org/copyleft/gpl.txt .
!!
!! INPUTS
!!  op=ctqmc
!!
!! OUTPUT
!!
!! SIDE EFFECTS
!!
!! NOTES
!!
!! PARENTS
!!  Will be filled automatically by the parent script
!!
!! CHILDREN
!!  Will be filled automatically by the parent script
!!
!! SOURCE

SUBROUTINE Ctqmc_allocateAll(op)

!Arguments ------------------------------------

!This section has been created automatically by the script Abilint (TD).
!Do not modify the following lines by hand.
#undef ABI_FUNC
#define ABI_FUNC 'Ctqmc_allocateAll'
!End of the abilint section

  TYPE(Ctqmc), INTENT(INOUT) :: op
!Local variables ------------------------------
  INTEGER                  :: flavors

  IF ( .NOT. op%para ) &
    CALL ERROR("Ctqmc_allocateAll : Ctqmc_setParameters never called  ")

  flavors = op%flavors

  DT_FREEIF(op%Greens)
  DT_MALLOC(op%Greens,(1:flavors))

  FREEIF(op%measN)
  MALLOC(op%measN,(1:4,1:flavors))
  op%measN = 0.d0

  FREEIF(op%measDE)
  MALLOC(op%measDE,(1:flavors,1:flavors) )
  op%measDE = 0.d0

  FREEIF(op%mu)
  MALLOC(op%mu,(1:flavors) )
  op%mu = 0.d0
END SUBROUTINE Ctqmc_allocateAll
!!***

!#ifdef CTCtqmc_ANALYSIS
!!****f* ABINIT/m_Ctqmc/Ctqmc_allocateOpt
!! NAME
!!  Ctqmc_allocateOpt
!!
!! FUNCTION
!!  allocate all option variables 
!!
!! COPYRIGHT
!!  Copyright (C) 2013-2014 ABINIT group (J. Bieder)
!!  This file is distributed under the terms of the
!!  GNU General Public License, see ~abinit/COPYING
!!  or http://www.gnu.org/copyleft/gpl.txt .
!!
!! INPUTS
!!  op=ctqmc
!!
!! OUTPUT
!!
!! SIDE EFFECTS
!!
!! NOTES
!!
!! PARENTS
!!  Will be filled automatically by the parent script
!!
!! CHILDREN
!!  Will be filled automatically by the parent script
!!
!! SOURCE

SUBROUTINE Ctqmc_allocateOpt(op)

!Arguments ------------------------------------

!This section has been created automatically by the script Abilint (TD).
!Do not modify the following lines by hand.
#undef ABI_FUNC
#define ABI_FUNC 'Ctqmc_allocateOpt'
!End of the abilint section

  TYPE(Ctqmc), INTENT(INOUT) :: op
!Local variables ------------------------------
  INTEGER :: i
  INTEGER :: j
  INTEGER :: k

  IF ( .NOT. op%para ) &
    CALL ERROR("Ctqmc_allocateOpt : Ctqmc_setParameters never called  ")

  IF ( op%opt_analysis .EQ. 1 ) THEN
    FREEIF(op%measCorrelation)
    MALLOC(op%measCorrelation,(1:op%samples+1,1:3,1:op%flavors))
    op%measCorrelation = 0.d0
  END IF

  IF ( op%opt_order .GT. 0 ) THEN
    FREEIF(op%measPerturbation)
    MALLOC(op%measPerturbation,(1:op%opt_order,1:op%flavors))
    op%measPerturbation = 0.d0
  END IF

  IF ( op%opt_noise .EQ. 1 ) THEN
    IF ( ALLOCATED(op%measNoiseG) ) THEN
      DO i=1,2
        DO j = 1, op%flavors
          DO k= 1, op%samples+1
            CALL Vector_destroy(op%measNoiseG(k,j,i))
          END DO
        END DO
      END DO
      DT_FREE(op%measNoiseG)
    END IF
    DT_MALLOC(op%measNoiseG,(1:op%samples+1,1:op%flavors,1:2))
    !DO i=1,2
      DO j = 1, op%flavors
        DO k= 1, op%samples+1
          CALL Vector_init(op%measNoiseG(k,j,1),CTQMC_SLICE1)
        END DO
      END DO
      DO j = 1, op%flavors
        DO k= 1, op%samples+1
          CALL Vector_init(op%measNoiseG(k,j,2),CTQMC_SLICE1*CTQMC_SLICE2+1) ! +1 pour etre remplacer ceil
        END DO
      END DO
    !END DO
    FREEIF(op%abNoiseG)
    MALLOC(op%aBNoiseG,(1:2,1:op%samples+1,op%flavors))
    op%abNoiseG = 0.d0
  END IF

  IF (op%opt_spectra .GE. 1 ) THEN
    FREEIF(op%density)
    !MALLOC(op%density,(1:op%thermalization,1:op%flavors))
    i = CEILING(DBLE(op%thermalization+op%sweeps)/DBLE(op%measurements*op%opt_spectra))
    MALLOC(op%density,(1:op%flavors+1,1:i))
    op%density = 0.d0
  END IF
!#endif
END SUBROUTINE Ctqmc_allocateOpt
!!***

SUBROUTINE Ctqmc_setG0wFile(op,istream,opt_fk)


!This section has been created automatically by the script Abilint (TD).
!Do not modify the following lines by hand.
#undef ABI_FUNC
#define ABI_FUNC 'Ctqmc_setG0wFile'
!End of the abilint section


#ifdef HAVE_MPI1
include 'mpif.h'
#endif
  TYPE(Ctqmc)                    , INTENT(INOUT) :: op
  INTEGER                        , INTENT(IN   ) :: istream
  INTEGER                        , INTENT(IN   ) :: opt_fk
  COMPLEX(KIND=8), DIMENSION(:,:), ALLOCATABLE   :: Gomega
  INTEGER                                        :: flavors
  INTEGER                                        :: iflavor
  INTEGER                                        :: iomega
#ifdef HAVE_MPI
  INTEGER                                        :: ierr
#endif

  IF ( .NOT. op%para ) &
    CALL ERROR("Ctqmc_setG0wFile : Ctqmc_setParameters never called   ") 

  flavors = op%flavors

  op%Wmax = op%samples
  MALLOC(Gomega,(1:op%Wmax,1:flavors))

  IF ( op%rank .EQ. 0 ) THEN
    DO iomega=1, op%Wmax
      READ(istream,*) (Gomega(iomega,iflavor),iflavor=1,flavors)
    END DO
  END IF

#ifdef HAVE_MPI
  CALL MPI_Bcast(Gomega, op%Wmax*flavors, MPI_DOUBLE_COMPLEX, 0,    &
                 op%MY_COMM, ierr)
#endif

  CALL Ctqmc_setG0wTab(op,Gomega,opt_fk)
  FREE(Gomega)

END SUBROUTINE Ctqmc_setG0wFile
!!***

!!****f* ABINIT/m_Ctqmc/Ctqmc_setG0wTab
!! NAME
!!  Ctqmc_setG0wTab
!!
!! FUNCTION
!!  Set Gow from input array
!!
!! COPYRIGHT
!!  Copyright (C) 2013-2014 ABINIT group (J. Bieder)
!!  This file is distributed under the terms of the
!!  GNU General Public License, see ~abinit/COPYING
!!  or http://www.gnu.org/copyleft/gpl.txt .
!!
!! INPUTS
!!  op=ctqmc
!!  Gomega=G0w
!!  opt_fk=F is already inversed with out iwn
!!
!! OUTPUT
!!
!! SIDE EFFECTS
!!
!! NOTES
!!
!! PARENTS
!!  Will be filled automatically by the parent script
!!
!! CHILDREN
!!  Will be filled automatically by the parent script
!!
!! SOURCE

SUBROUTINE Ctqmc_setG0wTab(op,Gomega,opt_fk)

!Arguments ------------------------------------

!This section has been created automatically by the script Abilint (TD).
!Do not modify the following lines by hand.
#undef ABI_FUNC
#define ABI_FUNC 'Ctqmc_setG0wTab'
!End of the abilint section

  TYPE(Ctqmc), INTENT(INOUT)                      :: op
  COMPLEX(KIND=8), DIMENSION(:,:), INTENT(IN ) :: Gomega
  INTEGER                         , INTENT(IN ) :: opt_fk
!Local variable -------------------------------
  DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE :: F

  IF ( .NOT. op%para ) &
    CALL ERROR("Ctqmc_setG0wTab : Ctqmc_setParameters never called    ") 

  MALLOC(F,(1:op%samples+1,1:op%flavors))
  CALL Ctqmc_computeF(op,Gomega, F, opt_fk)  ! mu is changed
  CALL BathOperator_setF(op%Bath, F)
  !CALL BathOperator_printF(op%Bath)
  FREE(F)

  op%inF = .TRUE.
  op%set = .TRUE. 

END SUBROUTINE Ctqmc_setG0wTab
!!***

!SUBROUTINE Ctqmc_setFwK(op,Gomega)
!  COMPLEX*16      , DIMENSION(:,:), INTENT(IN ) :: Gomega
!  DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE :: F
!
!  IF ( .NOT. op%para ) &
!    CALL ERROR("Ctqmc_setG0wTab : Ctqmc_setParameters never called    ") 
!
!  MALLOC(F,(1:op%samples+1,1:op%flavors))
!  CALL Ctqmc_computeFK(op,Gomega, op%Wmax, F)  ! mu is changed
!  CALL BathOperator_setF(op%Bath, F)
!  CALL BathOperator_printF(op%Bath)
!  FREE(F)
!
!  op%inF = .TRUE.
!  op%set = .TRUE. 
!
!END SUBROUTINE Ctqmc_setFwK
!!***

!SUBROUTINE Ctqmc_setBand(op, mu, U)
!  DOUBLE PRECISION, INTENT(IN   ) :: mu
!  DOUBLE PRECISION, INTENT(IN   ) :: U
!
!  IF ( .NOT. op%para ) &
!    CALL ERROR("Ctqmc_setBand : Ctqmc_setParameters never called      ")
!
!  CALL ImpurityOperator_setU(op%Impurity, U, 0.d0)
!  !op%mu = mu + op%Impurity%shift_mu
!END SUBROUTINE Ctqmc_setBand
!!***

!!****f* ABINIT/m_Ctqmc/Ctqmc_setU
!! NAME
!!  Ctqmc_setU
!!
!! FUNCTION
!!  set the interaction matrix
!!
!! COPYRIGHT
!!  Copyright (C) 2013-2014 ABINIT group (J. Bieder)
!!  This file is distributed under the terms of the
!!  GNU General Public License, see ~abinit/COPYING
!!  or http://www.gnu.org/copyleft/gpl.txt .
!!
!! INPUTS
!!  op=ctqmc
!!  matU=interaction matrix
!!
!! OUTPUT
!!
!! SIDE EFFECTS
!!
!! NOTES
!!
!! PARENTS
!!  Will be filled automatically by the parent script
!!
!! CHILDREN
!!  Will be filled automatically by the parent script
!!
!! SOURCE

SUBROUTINE Ctqmc_setU(op,matU)

!Arguments ------------------------------------

!This section has been created automatically by the script Abilint (TD).
!Do not modify the following lines by hand.
#undef ABI_FUNC
#define ABI_FUNC 'Ctqmc_setU'
!End of the abilint section

  TYPE(Ctqmc), INTENT(INOUT) ::op
!Local variables ------------------------------
  DOUBLE PRECISION, DIMENSION(:,:), INTENT(IN) :: matU

  IF ( SIZE(matU) .NE. op%flavors*op%flavors ) &
    CALL ERROR("Ctqmc_setU : Wrong interaction matrix (size)        ")

  CALL ImpurityOperator_setUmat(op%Impurity, matU)
  op%setU = .TRUE.
END SUBROUTINE Ctqmc_setU
!!***

!!****f* ABINIT/m_Ctqmc/Ctqmc_clear
!! NAME
!!  Ctqmc_clear
!!
!! FUNCTION
!!  clear a ctqmc run
!!
!! COPYRIGHT
!!  Copyright (C) 2013-2014 ABINIT group (J. Bieder)
!!  This file is distributed under the terms of the
!!  GNU General Public License, see ~abinit/COPYING
!!  or http://www.gnu.org/copyleft/gpl.txt .
!!
!! INPUTS
!!  op=ctqmc
!!
!! OUTPUT
!!
!! SIDE EFFECTS
!!
!! NOTES
!!
!! PARENTS
!!  Will be filled automatically by the parent script
!!
!! CHILDREN
!!  Will be filled automatically by the parent script
!!
!! SOURCE

SUBROUTINE Ctqmc_clear(op)

!Arguments ------------------------------------

!This section has been created automatically by the script Abilint (TD).
!Do not modify the following lines by hand.
#undef ABI_FUNC
#define ABI_FUNC 'Ctqmc_clear'
!End of the abilint section

  TYPE(Ctqmc), INTENT(INOUT) :: op
!Local variables ------------------------------
  INTEGER :: i
  INTEGER :: j
  INTEGER :: k

  op%measN(1,:) = 0.d0
  op%measN(2,:) = 0.d0
  !Do not set measN(3,:) to 0 to avoid erasing N between therm and ctqmc
  op%measN(4,:) = 0.d0
  op%measDE = 0.d0
!  op%seg_added    = 0.d0
!  op%anti_added   = 0.d0
!  op%seg_removed  = 0.d0
!  op%anti_removed = 0.d0
!  op%seg_sign     = 0.d0
!  op%anti_sign    = 0.d0
  op%stats(:)     = 0.d0
  op%swap         = 0.d0
  op%runTime      = 0.d0
  op%modGlobalMove(2) = 0 
  CALL Vector_clear(op%measNoise(1))
  CALL Vector_clear(op%measNoise(2))
!#ifdef CTCtqmc_CHECK
  op%errorImpurity = 0.d0
  op%errorBath     = 0.d0
!#endif
  DO j = 1, op%flavors
    CALL GreenHyb_clear(op%Greens(j))
  END DO
!#ifdef CTCtqmc_ANALYSIS
  IF ( op%opt_analysis .EQ. 1 .AND. ALLOCATED(op%measCorrelation) ) &    
    op%measCorrelation = 0.d0 
  IF ( op%opt_order .GT. 0 .AND. ALLOCATED(op%measPerturbation) ) &
    op%measPerturbation = 0.d0
  IF ( op%opt_noise .EQ. 1 .AND. ALLOCATED(op%measNoiseG) ) THEN
    DO i=1,2
      DO j = 1, op%flavors
        DO k= 1, op%samples+1
          CALL Vector_clear(op%measNoiseG(k,j,i))
        END DO
      END DO
    END DO
  END IF
!#endif
END SUBROUTINE Ctqmc_clear
!!***

!!****f* ABINIT/m_Ctqmc/Ctqmc_reset
!! NAME
!!  Ctqmc_reset
!!
!! FUNCTION
!!  reset a ctqmc simulation
!!
!! COPYRIGHT
!!  Copyright (C) 2013-2014 ABINIT group (J. Bieder)
!!  This file is distributed under the terms of the
!!  GNU General Public License, see ~abinit/COPYING
!!  or http://www.gnu.org/copyleft/gpl.txt .
!!
!! INPUTS
!!  op=ctqmc
!!
!! OUTPUT
!!
!! SIDE EFFECTS
!!
!! NOTES
!!
!! PARENTS
!!  Will be filled automatically by the parent script
!!
!! CHILDREN
!!  Will be filled automatically by the parent script
!!
!! SOURCE

SUBROUTINE Ctqmc_reset(op)

!Arguments ------------------------------------

!This section has been created automatically by the script Abilint (TD).
!Do not modify the following lines by hand.
#undef ABI_FUNC
#define ABI_FUNC 'Ctqmc_reset'
!End of the abilint section

  TYPE(Ctqmc), INTENT(INOUT) :: op
!Local variables ------------------------------
  INTEGER                  :: iflavor
  DOUBLE PRECISION         :: sweeps

  DO iflavor = 1, op%flavors
    CALL GreenHyb_reset(op%Greens(iflavor))
  END DO
  CALL Ctqmc_clear(op)
  CALL ImpurityOperator_reset(op%Impurity)
  CALL BathOperator_reset    (op%Bath)
  op%measN(3,:) = 0.d0
  !complete restart -> measN=0
  op%done = .FALSE.
  op%set  = .FALSE.
  op%inF  = .FALSE.
  op%opt_movie = 0
  op%opt_analysis = 0
  op%opt_order = 0
  op%opt_check = 0
  op%opt_noise = 0
  op%opt_spectra = 0
  op%opt_levels = 0
  sweeps = DBLE(op%sweeps)*DBLE(op%size)
  CALL Ctqmc_setSweeps(op, sweeps)
!#ifdef HAVE_MPI
!  CALL MPI_BARRIER(op%MY_COMM,iflavor)
!  IF ( op%rank .EQ. 0 ) &
!#endif
!  WRITE(op%ostream,'(A9)') "QMC reset"
!  CALL FLUSH(op%ostream)
END SUBROUTINE Ctqmc_reset
!!***

!!****f* ABINIT/m_Ctqmc/Ctqmc_setMu
!! NAME
!!  Ctqmc_setMu
!!
!! FUNCTION
!!  impose energy levels
!!
!! COPYRIGHT
!!  Copyright (C) 2013-2014 ABINIT group (J. Bieder)
!!  This file is distributed under the terms of the
!!  GNU General Public License, see ~abinit/COPYING
!!  or http://www.gnu.org/copyleft/gpl.txt .
!!
!! INPUTS
!!  op=ctqmc
!!  levels=energy levels vector
!!
!! OUTPUT
!!  argout(sizeout)=description
!!
!! SIDE EFFECTS
!!
!! NOTES
!!
!! PARENTS
!!  Will be filled automatically by the parent script
!!
!! CHILDREN
!!  Will be filled automatically by the parent script
!!
!! SOURCE

SUBROUTINE Ctqmc_setMu(op, levels)

!Arguments ------------------------------------

!This section has been created automatically by the script Abilint (TD).
!Do not modify the following lines by hand.
#undef ABI_FUNC
#define ABI_FUNC 'Ctqmc_setMu'
!End of the abilint section

  TYPE(Ctqmc)                   , INTENT(INOUT) :: op
  DOUBLE PRECISION, DIMENSION(:), INTENT(IN   ) :: levels

  IF ( op%flavors .NE. SIZE(levels,1) ) &
    CALL WARNALL("Ctqmc_setMu : Taking energy levels from weiss G(iw)")

  op%mu(:)=-levels  ! levels = \epsilon_j - \mu
  !op%mu =\tilde{\mu} = \mu -\epsilon_j
  op%opt_levels = 1
END SUBROUTINE Ctqmc_setMu
!!***

!!****f* ABINIT/m_Ctqmc/Ctqmc_computeF
!! NAME
!!  Ctqmc_computeF
!!
!! FUNCTION
!!  Compute the hybridization function
!!
!! COPYRIGHT
!!  Copyright (C) 2013-2014 ABINIT group (J. Bieder)
!!  This file is distributed under the terms of the
!!  GNU General Public License, see ~abinit/COPYING
!!  or http://www.gnu.org/copyleft/gpl.txt .
!!
!! INPUTS
!!  op=ctqmc
!!  Gomega=G0 to compute F
!!  opt_fk=What is Gomega
!!
!! OUTPUT
!!  F=hybridization function
!!
!! SIDE EFFECTS
!!
!! NOTES
!!
!! PARENTS
!!  Will be filled automatically by the parent script
!!
!! CHILDREN
!!  Will be filled automatically by the parent script
!!
!! SOURCE

SUBROUTINE Ctqmc_computeF(op, Gomega, F, opt_fk)

!Arguments ------------------------------------

!This section has been created automatically by the script Abilint (TD).
!Do not modify the following lines by hand.
#undef ABI_FUNC
#define ABI_FUNC 'Ctqmc_computeF'
!End of the abilint section

  TYPE(Ctqmc)                       , INTENT(INOUT) :: op
  COMPLEX(KIND=8), DIMENSION(:,:), INTENT(IN   ) :: Gomega
  !INTEGER                         , INTENT(IN   ) :: Wmax
  DOUBLE PRECISION, DIMENSION(:,:), INTENT(INOUT) :: F
  INTEGER                         , INTENT(IN   ) :: opt_fk
!Local variables ------------------------------
  INTEGER                                         :: flavors
  INTEGER                                         :: samples
  INTEGER                                         :: iflavor
  INTEGER                                         :: iomega
  INTEGER                                         :: itau
  DOUBLE PRECISION                                :: pi_invBeta
  DOUBLE PRECISION                                :: K
  DOUBLE PRECISION                                :: re
  DOUBLE PRECISION                                :: im
  COMPLEX(KIND=8), DIMENSION(:,:), ALLOCATABLE   :: F_omega
  TYPE(GreenHyb)                                     :: F_tmp

  flavors    = op%flavors

  samples    = op%samples
  pi_invBeta = ACOS(-1.d0) / op%beta
  op%Wmax=SIZE(Gomega,1)

  IF ( op%have_MPI .EQV. .TRUE. ) THEN
    CALL GreenHyb_init(F_tmp,samples,op%beta,MY_COMM=op%MY_COMM)
  ELSE
    CALL GreenHyb_init(F_tmp,samples,op%beta)
  END IF
!  K = op%mu

  MALLOC(F_omega,(1:op%Wmax,1:flavors))

  !IF ( op%rank .EQ. 0 ) &
    !OPEN(UNIT=9876,FILE="K.dat",POSITION="APPEND")
  IF ( opt_fk .EQ. 0 ) THEN
    DO iflavor = 1, flavors
      DO iomega=1,op%Wmax
        re = REAL(Gomega(iomega,iflavor))
        im = AIMAG(Gomega(iomega,iflavor))
        F_omega(iomega,iflavor) = CMPLX(-re/(re*re+im*im),im/(re*re+im*im),8)
      END DO
    END DO
    !F_omega = CMPLX(-1.d0,0,8)/Gomega
  ELSE
    F_omega = Gomega
  END IF

  DO iflavor = 1, flavors
    IF ( op%opt_levels .EQ. 1 ) THEN
      K = op%mu(iflavor)
    ELSE
      K = -REAL(F_omega(op%Wmax, iflavor))
!    op%mu = K
      op%mu(iflavor) = K 
    END IF
    !IF ( op%rank .EQ. 0 ) &
    !WRITE(9876,'(I4,2E22.14)') iflavor, K, REAL(-F_omega(op%Wmax, iflavor))
    !IF(op%rank .EQ.0) &
    !WRITE(op%ostream,*) "CTQMC K, op%mu = ",K,op%mu(iflavor)
    !WRITE(op%ostream,*) "CTQMC beta     = ",op%beta
    IF ( opt_fk .EQ. 0 ) THEN
      DO iomega = 1, op%Wmax
        re = REAL(F_omega(iomega,iflavor))
        im = AIMAG(F_omega(iomega,iflavor))
        F_omega(iomega,iflavor) = CMPLX(re + K, im + (2.d0*DBLE(iomega)-1.d0) * pi_invBeta, 8)
        !if(iflavor==1.and.op%rank==0) then
          !write(224,*) (2.d0*DBLE(iomega)-1.d0) * pi_invBeta, real(F_omega(iomega,iflavor)),imag(F_omega(iomega,iflavor))
          !write(225,*) (2.d0*DBLE(iomega)-1.d0) * pi_invBeta, real(Gomega(iomega, iflavor)),imag(Gomega(iomega, iflavor))
        !end if 
      END DO
    ELSE
      DO iomega = 1, op%Wmax
        F_omega(iomega,iflavor) = F_omega(iomega,iflavor) &
                    + CMPLX(K, 0.d0, 8)
        !if(iflavor==1.and.op%rank==0) then
          !write(224,*) (2.d0*DBLE(iomega)-1.d0) * pi_invBeta, real(F_omega(iomega,iflavor)),imag(F_omega(iomega,iflavor))
          !write(225,*) (2.d0*DBLE(iomega)-1.d0) * pi_invBeta, real(Gomega(iomega, iflavor)),imag(Gomega(iomega, iflavor))
        !end if 
      END DO
    END IF
    K = REAL(CMPLX(0,(2.d0*DBLE(op%Wmax)-1.d0)*pi_invBeta,8)*F_omega(op%Wmax,iflavor))
    CALL GreenHyb_setMuD1(op%Greens(iflavor),op%mu(iflavor),K)
    CALL GreenHyb_setOperW(F_tmp,F_omega(:,iflavor))
    !CALL GreenHyb_backFourier(F_tmp,F_omega(:,iflavor))
    CALL GreenHyb_backFourier(F_tmp)
    F(1:samples+1,iflavor) = (/ (-F_tmp%oper(samples+1-itau),itau=0,samples) /)
  END DO
  FREE(F_omega)
  CALL GreenHyb_destroy(F_tmp)
END SUBROUTINE Ctqmc_computeF
!!***

!SUBROUTINE Ctqmc_computeFK(op, Gomega, Wmax, F)
!  COMPLEX*16      , DIMENSION(:,:), INTENT(IN   ) :: Gomega
!  INTEGER                         , INTENT(IN   ) :: Wmax
!  DOUBLE PRECISION, DIMENSION(:,:), INTENT(INOUT) :: F
!  INTEGER                                         :: flavors
!  INTEGER                                         :: samples
!  INTEGER                                         :: iflavor
!  INTEGER                                         :: iomega
!  INTEGER                                         :: itau
!  DOUBLE PRECISION                                :: pi_invBeta
!  DOUBLE PRECISION                                :: K
!  COMPLEX*16      , DIMENSION(:,:), ALLOCATABLE   :: F_omega
!  TYPE(GreenHyb)                                     :: F_tmp
!
!  flavors    = op%flavors
!
!  samples    = op%samples
!  pi_invBeta = ACOS(-1.d0) / op%beta
!
!  IF ( op%have_MPI .EQV. .TRUE. ) THEN
!    CALL GreenHyb_init(F_tmp,samples,op%beta,op%MY_COMM)
!  ELSE
!    CALL GreenHyb_init(F_tmp,samples,op%beta)
!  END IF
!!  K = op%mu
!
!  MALLOC(F_omega,(1:Wmax,1:flavors))
!
!  DO iflavor = 1, flavors
!    K = REAL(Gomega(Wmax, iflavor))
!    WRITE(op%ostream,*) "CTQMC K, op%mu = ",K,op%mu
!    WRITE(op%ostream,*) "CTQMC beta     = ",op%beta
!    op%mu(iflavor) = K 
!    DO iomega = 1, Wmax
!      F_omega(iomega,iflavor) = Gomega(iomega,iflavor) &
!                  - CMPLX(K, 0.d0, 8)
!      !if(iflavor==1.and.op%rank==0) then
!        !write(224,*) (2.d0*DBLE(iomega)-1.d0) * pi_invBeta, real(F_omega(iomega,iflavor)),imag(F_omega(iomega,iflavor))
!        !write(225,*) (2.d0*DBLE(iomega)-1.d0) * pi_invBeta, real(Gomega(iomega, iflavor)),imag(Gomega(iomega, iflavor))
!      !end if 
!    END DO
!    CALL GreenHyb_backFourier(F_tmp,F_omega(:,iflavor))
!    F(1:samples+1,iflavor) = (/ (-F_tmp%oper(samples+1-itau),itau=0,samples) /)
!  END DO
!  FREE(F_omega)
!  CALL GreenHyb_destroy(F_tmp)
!END SUBROUTINE Ctqmc_computeFK
!!***

!!****f* ABINIT/m_Ctqmc/Ctqmc_run
!! NAME
!!  Ctqmc_run
!!
!! FUNCTION
!!  set all options and run a simulation
!!
!! COPYRIGHT
!!  Copyright (C) 2013-2014 ABINIT group (J. Bieder)
!!  This file is distributed under the terms of the
!!  GNU General Public License, see ~abinit/COPYING
!!  or http://www.gnu.org/copyleft/gpl.txt .
!!
!! INPUTS
!!  op=ctqmc
!!  opt_order=maximal perturbation order to scope
!!  opt_movie=draw a movie of the simulation
!!  opt_analysis=compute correlation functions
!!  opt_check=check fast calculations
!!  opt_noise=compute noise for green function
!!  opt_spectra=fourier transform of the time evolution of the number of electrons
!!  opt_gMove=steps without global move
!!
!! OUTPUT
!!
!! SIDE EFFECTS
!!
!! NOTES
!!
!! PARENTS
!!  Will be filled automatically by the parent script
!!
!! CHILDREN
!!  Will be filled automatically by the parent script
!!
!! SOURCE

SUBROUTINE Ctqmc_run(op,opt_order,opt_movie,opt_analysis,opt_check,opt_noise,opt_spectra,opt_gMove)


!This section has been created automatically by the script Abilint (TD).
!Do not modify the following lines by hand.
#undef ABI_FUNC
#define ABI_FUNC 'Ctqmc_run'
!End of the abilint section


#ifdef HAVE_MPI1
include 'mpif.h'
#endif
!Arguments ------------------------------------
  TYPE(Ctqmc), INTENT(INOUT)           :: op
  INTEGER, OPTIONAL, INTENT(IN   )  :: opt_order
  INTEGER, OPTIONAL, INTENT(IN   )  :: opt_movie
  INTEGER, OPTIONAL, INTENT(IN   )  :: opt_analysis
  INTEGER, OPTIONAL, INTENT(IN   )  :: opt_check
  INTEGER, OPTIONAL, INTENT(IN   )  :: opt_noise
  INTEGER, OPTIONAL, INTENT(IN   )  :: opt_spectra
  INTEGER, OPTIONAL, INTENT(IN   )  :: opt_gMove
!Local variables ------------------------------
#ifdef HAVE_MPI
  INTEGER                            :: ierr
#endif
!#ifdef CTCtqmc_MOVIE
  INTEGER                            :: ilatex
  CHARACTER(LEN=4)                   :: Cchar
!#endif
  DOUBLE PRECISION                   :: estimatedTime

  IF ( .NOT. op%set  ) &
    CALL ERROR("Ctqmc_run : QMC not set up                          ")
  IF ( .NOT. op%setU ) &
    CALL ERROR("Ctqmc_run : QMC does not have a U matrix            ")


! OPTIONS of the run
  IF ( PRESENT( opt_check ) ) THEN
    op%opt_check = opt_check
    CALL ImpurityOperator_doCheck(op%Impurity,opt_check)
    CALL BathOperator_doCheck(op%Bath,opt_check)
  END IF
  IF ( PRESENT( opt_movie ) ) &
    op%opt_movie = opt_movie
  IF ( PRESENT( opt_analysis ) ) &
    op%opt_analysis = opt_analysis
  IF ( PRESENT ( opt_order ) ) &
    op%opt_order = opt_order 
  IF ( PRESENT ( opt_noise ) ) THEN
    op%opt_noise = opt_noise 
  END IF
  IF ( PRESENT ( opt_spectra ) ) &
    op%opt_spectra = opt_spectra

  op%modGlobalMove(1) = op%sweeps+1 ! No Global Move
  op%modGlobalMove(2) = 0
  IF ( PRESENT ( opt_gMove ) ) THEN
    IF ( opt_gMove .LE. 0 .OR. opt_gMove .GT. op%sweeps ) THEN
      op%modGlobalMove(1) = op%sweeps+1
      CALL WARNALL("Ctqmc_run : global moves option is <= 0 or > sweeps/cpu -> No global Moves")
    ELSE 
      op%modGlobalMove(1) = opt_gMove 
    END IF
  END IF

  CALL Ctqmc_allocateOpt(op)
  
!#ifdef CTCtqmc_MOVIE  
  ilatex = 0
  IF ( op%opt_movie .EQ. 1 ) THEN
    Cchar ="0000"
    WRITE(Cchar,'(I4)') op%rank 
    ilatex = 87+op%rank
    OPEN(UNIT=ilatex, FILE="Movie_"//TRIM(ADJUSTL(Cchar))//".tex")
    WRITE(ilatex,'(A)') "\documentclass{beamer}"
    WRITE(ilatex,'(A)') "\usepackage{color}"
    WRITE(ilatex,'(A)') "\setbeamersize{sidebar width left=0pt}"
    WRITE(ilatex,'(A)') "\setbeamersize{sidebar width right=0pt}"
    WRITE(ilatex,'(A)') "\setbeamersize{text width left=0pt}"
    WRITE(ilatex,'(A)') "\setbeamersize{text width right=0pt}"
    WRITE(ilatex,*) 
    WRITE(ilatex,'(A)') "\begin{document}"
    WRITE(ilatex,*) 
  END IF
!#endif

  IF ( op%rank .EQ. 0 ) THEN
    WRITE(op%ostream,'(A29)') "Starting QMC (Thermalization)"
  END IF
  
  !=================================
  ! STARTING THERMALIZATION 
  !=================================
  CALL Ctqmc_loop(op,op%thermalization,ilatex)
  !=================================
  ! ENDING   THERMALIZATION 
  !=================================

  estimatedTime = op%runTime
#ifdef HAVE_MPI
  CALL MPI_REDUCE(op%runTime, estimatedTime, 1, MPI_DOUBLE_PRECISION, MPI_MAX, &
             0, op%MY_COMM, ierr)
#endif

  IF ( op%rank .EQ. 0 ) THEN
    WRITE(op%ostream,'(A26,I6,A11)') "Thermalization done in    ", CEILING(estimatedTime), "    seconds"
    WRITE(op%ostream,'(A25,I7,A15,I5,A5)') "The QMC should run in    ", &
           CEILING(estimatedTime*DBLE(op%sweeps)/DBLE(op%thermalization)),&
                        "    seconds on ", op%size, " CPUs"
  END IF

  !=================================
  ! CLEANING CTQMC          
  !=================================
  CALL Ctqmc_clear(op)

  !=================================
  ! STARTING CTQMC          
  !=================================
  CALL Ctqmc_loop(op,op%sweeps,ilatex)
  !=================================
  ! ENDING   CTQMC          
  !=================================

  IF ( op%opt_movie .EQ. 1 ) THEN
    WRITE(ilatex,*) ""
    WRITE(ilatex,'(A14)') "\end{document}"
    CLOSE(ilatex)
  END IF

  op%done     = .TRUE.

END SUBROUTINE Ctqmc_run
!!***

!!****f* ABINIT/m_Ctqmc/Ctqmc_loop
!! NAME
!!  Ctqmc_loop
!!
!! FUNCTION
!!  Definition the main loop of the CT-QMC
!!
!! COPYRIGHT
!!  Copyright (C) 2013 ABINIT group (J. Bieder)
!!  This file is distributed under the terms of the
!!  GNU General Public License, see ~abinit/COPYING
!!  or http://www.gnu.org/copyleft/gpl.txt .
!!
!! INPUTS
!!  op=ctqmc
!!  itotal=number of sweeps to perform : thermalization or sweeps
!!  ilatex=unit of file to write movie if so
!!
!! OUTPUT
!!
!! SIDE EFFECTS
!!
!! NOTES
!!
!! PARENTS
!!  Will be filled automatically by the parent script
!!
!! CHILDREN
!!  Will be filled automatically by the parent script
!!
!! SOURCE

SUBROUTINE Ctqmc_loop(op,itotal,ilatex)

!Arguments ------------------------------------

!This section has been created automatically by the script Abilint (TD).
!Do not modify the following lines by hand.
#undef ABI_FUNC
#define ABI_FUNC 'Ctqmc_loop'
!End of the abilint section

  TYPE(Ctqmc), INTENT(INOUT)         :: op
  INTEGER    , INTENT(IN   )         :: itotal
  INTEGER    , INTENT(IN   )         :: ilatex
!Local variables ------------------------------
  LOGICAL                            :: updated 
  LOGICAL                            :: updated_seg
  LOGICAL, DIMENSION(:), ALLOCATABLE :: updated_swap

  INTEGER                            :: flavors
  INTEGER                            :: measurements
  INTEGER                            :: modNoise1
  INTEGER                            :: modNoise2
  INTEGER                            :: modGlobalMove
  INTEGER                            :: sp1
  INTEGER                            :: itau   
  INTEGER                            :: ind
  INTEGER                            :: endDensity
  INTEGER                            :: indDensity
  INTEGER                            :: swapUpdate1
  INTEGER                            :: swapUpdate2
  INTEGER                            :: old_percent
  INTEGER                            :: new_percent
  INTEGER                            :: ipercent
  INTEGER                            :: iflavor
  INTEGER                            :: isweep

  DOUBLE PRECISION                   :: cpu_time1
  DOUBLE PRECISION                   :: cpu_time2
  DOUBLE PRECISION                   :: NRJ_old1
  DOUBLE PRECISION                   :: NRJ_old2
  DOUBLE PRECISION                   :: NRJ_new
  DOUBLE PRECISION                   :: total
  DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:) :: gtmp_old1
  DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:) :: gtmp_old2
  DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:) :: gtmp_new

  CALL CPU_TIME(cpu_time1)

  flavors        = op%flavors
  measurements   = op%measurements
  modNoise1      = op%modNoise1
  modNoise2      = op%modNoise2
  modGlobalMove  = op%modGlobalMove(1)
  sp1            = op%samples+1

  old_percent    = 0

  MALLOC(updated_swap,(1:flavors))
  updated_swap(:) = .FALSE.

  NRJ_old1  = 0.d0
  NRJ_old2  = 0.d0
  NRJ_new   = 0.d0

  MALLOC(gtmp_new,(1,1))
  gtmp_new  = 0.d0
  MALLOC(gtmp_old1,(1,1))
  gtmp_old1 = 0.d0
  MALLOC(gtmp_old2,(1,1))
  gtmp_old2 = 0.d0

  endDensity = SIZE(op%density,2)

  IF ( op%opt_noise .GT. 0 ) THEN
    FREEIF(gtmp_new)
    MALLOC(gtmp_new,(1:sp1,1:flavors))
    FREEIF(gtmp_old1)
    MALLOC(gtmp_old1,(1:sp1,1:flavors))
    FREEIF(gtmp_old2)
    MALLOC(gtmp_old2,(1:sp1,1:flavors))
  END IF

  IF ( op%rank .EQ. 0 ) THEN
    WRITE(op%ostream, '(1x,103A)') &
    "|----------------------------------------------------------------------------------------------------|"
    WRITE(op%ostream,'(1x,A)', ADVANCE="NO") "|"
  END IF

  total = DBLE(itotal)

  indDensity = 1
  DO isweep = 1, itotal
    DO iflavor = 1, flavors
      ImpurityOperator_QuickActivation(op%Impurity,iflavor)
      BathOperator_QuickActivation(op%Bath,iflavor)

      CALL Ctqmc_tryAddRemove(op,updated_seg)
      updated = updated_seg .OR.  updated_swap(iflavor)
      updated_swap(iflavor) = .FALSE.

      CALL GreenHyb_measHybrid(op%Greens(iflavor), op%Bath%M(iflavor), op%Impurity%Particles(iflavor), updated)
      CALL Ctqmc_measN        (op, iflavor, updated)
      IF ( op%opt_analysis .EQ. 1 ) &
        CALL Ctqmc_measCorrelation (op, iflavor)
      IF ( op%opt_order .GT. 0 ) &
        CALL Ctqmc_measPerturbation(op, iflavor)
    END DO

    IF ( MOD(isweep,modGlobalMove) .EQ. 0 ) THEN
      CALL Ctqmc_trySwap(op,swapUpdate1, swapUpdate2)
      IF ( swapUpdate1 .NE. 0 .AND. swapUpdate2 .NE. 0 ) THEN
        updated_swap(swapUpdate1) = .TRUE.
        updated_swap(swapUpdate2) = .TRUE.
      END IF
    END IF
    
    IF ( MOD(isweep,measurements) .EQ. 0 ) THEN
      CALL ImpurityOperator_measDE(op%Impurity,op%measDE)
      IF ( op%opt_spectra .GE. 1 .AND. MOD(isweep,measurements*op%opt_spectra) .EQ. 0 ) THEN
        op%density(1:flavors,indDensity) = op%measN(3,1:flavors)
        indDensity = indDensity+1
      END IF
    END IF

    IF ( MOD(isweep, modNoise1) .EQ. 0 ) THEN
      !modNext = isweep + modNoise2
      NRJ_new = (SUM(op%measDE(:,:))-op%measDE(1,1))*0.5d0 ! double occupation, avoid stat with 0 for U=J=0
      CALL Vector_pushBack(op%measNoise(1),NRJ_new - NRJ_old1)
      NRJ_old1 = NRJ_new

      !! Try to limit accumulation error
      CALL ImpurityOperator_cleanOverlaps(op%Impurity)

      IF ( op%opt_noise .EQ. 1 ) THEN
        DO iflavor = 1, flavors
          DO ind = 1, op%Greens(iflavor)%map%tail
            itau = op%Greens(iflavor)%map%listINT(ind)
            gtmp_new(itau,iflavor) = op%Greens(iflavor)%oper(itau) & 
                        +op%Greens(iflavor)%map%listDBLE(ind)*DBLE(op%Greens(iflavor)%factor)
          END DO
          DO itau = 1, sp1
            CALL Vector_pushBack(op%measNoiseG(itau,iflavor,1), gtmp_new(itau,iflavor) - gtmp_old1(itau,iflavor))
            gtmp_old1(itau,iflavor) = gtmp_new(itau,iflavor)
          END DO
        END DO
      END IF
    END IF

    IF ( MOD(isweep,modNoise2) .EQ. 0 ) THEN
      NRJ_new = (SUM(op%measDE(:,:))-op%measDE(1,1))*0.5d0 ! double occupation, avoid stat with 0 for U=J=0
      CALL Vector_pushBack(op%measNoise(2),NRJ_new - NRJ_old2)
      NRJ_old2 = NRJ_new
      IF ( op%opt_noise .EQ. 1 ) THEN
        DO iflavor = 1, flavors
          DO ind = 1, op%Greens(iflavor)%map%tail
            itau = op%Greens(iflavor)%map%listINT(ind)
            gtmp_new(itau,iflavor) = op%Greens(iflavor)%oper(itau) & 
                        +op%Greens(iflavor)%map%listDBLE(ind)*op%Greens(iflavor)%factor
          END DO
          DO itau = 1, sp1
            CALL Vector_pushBack(op%measNoiseG(itau,iflavor,2), gtmp_new(itau,iflavor) - gtmp_old2(itau,iflavor))
            gtmp_old2(itau,iflavor) = gtmp_new(itau,iflavor)
          END DO
        END DO 
      END IF

      IF ( op%rank .EQ. 0 ) THEN 
        new_percent = CEILING(DBLE(isweep)*100.d0/DBLE(itotal))
        DO ipercent = old_percent+1, new_percent 
          WRITE(op%ostream,'(A)',ADVANCE="NO") "-"
        END DO
        old_percent = new_percent
      END IF
    END IF

    IF ( op%opt_movie .EQ. 1 ) THEN
      WRITE(ilatex,'(A11,I9)') "%iteration ", isweep
      CALL ImpurityOperator_printLatex(op%Impurity,ilatex,isweep)
    END IF

  END DO

  IF ( op%rank .EQ. 0 ) THEN
    DO ipercent = old_percent+1, 100
      WRITE(op%ostream,'(A)',ADVANCE="NO") "-"
    END DO
    WRITE(op%ostream,'(A)') "|"
  END IF
 
  FREE(gtmp_new)
  FREE(gtmp_old1)
  FREE(gtmp_old2)
  FREE(updated_swap)

  IF ( op%opt_spectra .GE. 1 .AND. itotal .EQ. op%sweeps ) THEN
    IF ( endDensity .NE. indDensity-1 ) THEN
      op%density(:,endDensity) = -1.d0
    END IF
  END IF

  CALL CPU_TIME(cpu_time2)

  op%runTime = (cpu_time2 - cpu_time1)*1.05d0 ! facteur arbitraire de correction
END SUBROUTINE Ctqmc_loop
!!***

!!****f* ABINIT/m_Ctqmc/Ctqmc_tryAddRemove
!! NAME
!!  Ctqmc_tryAddRemove
!!
!! FUNCTION
!!  Try to add or remove a segment and an anti-segment
!!
!! COPYRIGHT
!!  Copyright (C) 2013-2014 ABINIT group (J. Bieder)
!!  This file is distributed under the terms of the
!!  GNU General Public License, see ~abinit/COPYING
!!  or http://www.gnu.org/copyleft/gpl.txt .
!!
!! INPUTS
!!  op=ctqmc
!!
!! OUTPUT
!!  updated=something changed
!!
!! SIDE EFFECTS
!!
!! NOTES
!!
!! PARENTS
!!  Will be filled automatically by the parent script
!!
!! CHILDREN
!!  Will be filled automatically by the parent script
!!
!! SOURCE

SUBROUTINE Ctqmc_tryAddRemove(op,updated)

!Arguments ------------------------------------

!This section has been created automatically by the script Abilint (TD).
!Do not modify the following lines by hand.
#undef ABI_FUNC
#define ABI_FUNC 'Ctqmc_tryAddRemove'
!End of the abilint section

  TYPE(Ctqmc)             , INTENT(INOUT) :: op
!  TYPE(BathOperator)    , INTENT(INOUT) :: Bath 
!  TYPE(ImpurityOperator), INTENT(INOUT) :: Impurity 
  LOGICAL               , INTENT(  OUT) :: updated
!Local variables ------------------------------
  INTEGER                               :: position
  INTEGER         , DIMENSION(1:2)     :: nature ! -2 for antiseg and 1 for seg
  INTEGER                               :: i! -2 for antiseg and 1 for seg
  DOUBLE PRECISION                      :: action
  DOUBLE PRECISION                      :: beta
  DOUBLE PRECISION                      :: time1
  DOUBLE PRECISION                      :: time2
  DOUBLE PRECISION                      :: time_avail
  DOUBLE PRECISION                      :: det_ratio
  DOUBLE PRECISION                      :: overlap
  DOUBLE PRECISION                      :: length
  DOUBLE PRECISION                      :: signe
  DOUBLE PRECISION                      :: tail
  DOUBLE PRECISION, DIMENSION(1:2)      :: CdagC_1

  IF ( .NOT. op%set ) &
    CALL ERROR("Ctqmc_trySegment : QMC not set                       ")

  nature(1) = CTQMC_SEGME
  nature(2) = CTQMC_ANTIS
  beta      = op%beta

  updated = .FALSE.
  tail  = DBLE(op%Impurity%particles(op%Impurity%activeFlavor)%tail)


  DO i = 1, 2
    signe = SIGN(1.d0,DBLE(nature(i))) 

    !CALL RANDOM_NUMBER(action)
    CALL OurRng(op%seed,action)

    IF ( action .LT. .5d0 ) THEN ! Ajout de segment
      !CALL RANDOM_NUMBER(time1)
      CALL OurRng(op%seed,time1)
      time1 = time1 * beta
      time_avail = ImpurityOperator_getAvailableTime(op%Impurity,time1,position) * signe
      IF ( time_avail .GT. 0.d0 ) THEN
        !CALL RANDOM_NUMBER(time2)
        CALL OurRng(op%seed,time2)
        IF ( time2 .EQ. 0.d0 ) THEN
          CALL OurRng(op%seed,time2) ! Prevent null segment
        END IF
        time2     = time1 + time2 * time_avail
!        CALL CdagC_init(CdagC_1,time1,time2)
          CdagC_1(Cdag_) = ((1.d0+signe)*time1+(1.d0-signe)*time2)*0.5d0
          CdagC_1(C_   ) = ((1.d0+signe)*time2+(1.d0-signe)*time1)*0.5d0
!        length    = CdagC_length(CdagC_1)
        length    = CdagC_1(C_   ) - CdagC_1(Cdag_)
        det_ratio = BathOperator_getDetAdd(op%Bath,CdagC_1,position,op%Impurity%particles(op%Impurity%activeFlavor))
        overlap   = ImpurityOperator_getNewOverlap(op%Impurity,CdagC_1)
        !CALL RANDOM_NUMBER(time1)
        CALL OurRng(op%seed,time1)
        IF ( det_ratio .LT. 0.d0 ) THEN
          det_ratio   = -det_ratio
          op%stats(nature(i)+CTQMC_DETSI) = op%stats(nature(i)+CTQMC_DETSI) + 1.d0
        END IF
        IF ( (time1 * (tail + 1.d0 )) &
             .LT. (beta * time_avail * det_ratio * DEXP(op%mu(op%Impurity%activeFlavor)*length + overlap) ) ) THEN
!          write(*,*) "before"
!          CALL ListCdagC_print(op%Impurity%particles(op%Impurity%activeFlavor),6)
          CALL ImpurityOperator_add(op%Impurity,CdagC_1,position)
!          write(*,*) "after "
!          CALL ListCdagC_print(op%Impurity%particles(op%Impurity%activeFlavor),6)
          CALL BathOperator_setMAdd(op%bath,op%Impurity%particles(op%Impurity%activeFlavor))
          op%stats(nature(i)+CTQMC_ADDED) = op%stats(nature(i)+CTQMC_ADDED)  + 1.d0
          updated = .TRUE. .OR. updated
          tail = tail + 1.d0
!          read(*,*) time1
        END IF 
      END IF 

    ELSE ! Supprimer un segment
      IF ( tail .GT. 0.d0 ) THEN
        !CALL RANDOM_NUMBER(time1)
        CALL OurRng(op%seed,time1)
        position = INT(((time1 * tail) + 1.d0) * signe )
        time_avail = ImpurityOperator_getAvailedTime(op%Impurity,position)
        det_ratio  = BathOperator_getDetRemove(op%Bath,position)
        CdagC_1    = ImpurityOperator_getSegment(op%Impurity,position)
!        length     = CdagC_length(CdagC_1)
        length     = CdagC_1(C_) - CdagC_1(Cdag_)
        overlap    = ImpurityOperator_getNewOverlap(op%Impurity,CdagC_1)
        !CALL RANDOM_NUMBER(time1)
        CALL OurRng(op%seed,time1)
        IF ( det_ratio .LT. 0.d0 ) THEN
          det_ratio   = -det_ratio
!          op%seg_sign = op%seg_sign + 1.d0
          op%stats(nature(i)+CTQMC_DETSI) = op%stats(nature(i)+CTQMC_DETSI) + 1.d0
        END IF
        IF ( (time1 * beta * time_avail * DEXP(op%mu(op%Impurity%activeFlavor)*length+overlap)) &
             .LT. (tail * det_ratio ) ) THEN
          CALL ImpurityOperator_remove(op%Impurity,position)
          CALL BathOperator_setMRemove(op%Bath,op%Impurity%particles(op%Impurity%activeFlavor))
          !op%seg_removed = op%seg_removed  + 1.d0
          op%stats(nature(i)+CTQMC_REMOV) = op%stats(nature(i)+CTQMC_REMOV)  + 1.d0
          updated = .TRUE. .OR. updated
          tail = tail -1.d0
        END IF
      END IF
    END IF
  END DO
END SUBROUTINE Ctqmc_tryAddRemove
!!***

!!SUBROUTINE Ctqmc_trySegment(op,updated)
!!!  TYPE(BathOperator)    , INTENT(INOUT) :: Bath 
!!!  TYPE(ImpurityOperator), INTENT(INOUT) :: Impurity 
!!  LOGICAL               , INTENT(INOUT) :: updated
!!  INTEGER                               :: position
!!  DOUBLE PRECISION                      :: action
!!  DOUBLE PRECISION                      :: beta
!!  DOUBLE PRECISION                      :: time1
!!  DOUBLE PRECISION                      :: time2
!!  DOUBLE PRECISION                      :: time_avail
!!  DOUBLE PRECISION                      :: det_ratio
!!  DOUBLE PRECISION                      :: overlap
!!  DOUBLE PRECISION                      :: length
!!  DOUBLE PRECISION                      :: tail
!!  DOUBLE PRECISION, DIMENSION(1:2)      :: CdagC_1
!!
!!  IF ( .NOT. op%set ) &
!!    CALL ERROR("Ctqmc_trySegment : QMC not set                       ")
!!
!!  beta     =  op%beta
!!  tail     = DBLE(op%Impurity%particles(op%Impurity%activeFlavor)%tail)
!!
!!  CALL RANDOM_NUMBER(action)
!!  
!!  updated = .FALSE.
!!
!!  IF ( action .LT. .5d0 ) THEN ! Ajout de segment
!!    CALL RANDOM_NUMBER(time1)
!!    time1 = time1 * beta
!!    time_avail = ImpurityOperator_getAvailableTime(op%Impurity,time1,position)
!!    IF ( time_avail .GT. 0.d0 ) THEN
!!      CALL RANDOM_NUMBER(time2)
!!      time2     = time1 + time2 * time_avail
!!!      CALL CdagC_init(CdagC_1,time1,time2)
!!      CdagC_1(Cdag_) = time1
!!      CdagC_1(C_   ) = time2
!!!      length    = CdagC_length(CdagC_1)
!!      length    = time2 - time1 
!!      det_ratio = BathOperator_getDetAdd(op%Bath,CdagC_1,position,op%Impurity%particles(op%Impurity%activeFlavor))
!!      overlap   = ImpurityOperator_getNewOverlap(op%Impurity,CdagC_1)
!!      CALL RANDOM_NUMBER(time1)
!!      IF ( det_ratio .LT. 0.d0 ) THEN
!!        det_ratio   = -det_ratio
!!        op%stats(CTQMC_SEGME+CTQMC_DETSI) = op%stats(CTQMC_SEGME+CTQMC_DETSI) + SIGN(1.d0,det_ratio)
!!      END IF
!!      IF ( (time1 * (tail + 1.d0 )) &
!!           .LT. (beta * time_avail * det_ratio * DEXP(op%mu(op%Impurity%activeFlavor)*length + overlap) ) ) THEN
!!        write(*,*) position
!!        CALL ImpurityOperator_add(op%Impurity,CdagC_1,position)
!!        CALL BathOperator_setMAdd(op%bath,op%Impurity%particles(op%Impurity%activeFlavor))
!!        op%stats(CTQMC_SEGME+CTQMC_ADDED) = op%stats(CTQMC_SEGME+CTQMC_ADDED)  + 1.d0
!!        updated = .TRUE.
!!      END IF 
!!    END IF 
!!
!!  ELSE ! Supprimer un segment
!!    IF ( tail .GT. 0.d0 ) THEN
!!      CALL RANDOM_NUMBER(time1)
!!      position = INT(time1 * tail) + 1
!!      time_avail = ImpurityOperator_getAvailedTime(op%Impurity,position)
!!      det_ratio  = BathOperator_getDetRemove(op%Bath,position)
!!      CdagC_1    = ImpurityOperator_getSegment(op%Impurity,position)
!!!      length     = CdagC_length(CdagC_1)
!!      length     = CdagC_1(C_) - CdagC_1(Cdag_)
!!      overlap    = ImpurityOperator_getNewOverlap(op%Impurity,CdagC_1)
!!      CALL RANDOM_NUMBER(time1)
!!      IF ( det_ratio .LT. 0.d0 ) THEN
!!        det_ratio   = -det_ratio
!!!        op%seg_sign = op%seg_sign + 1.d0
!!        op%stats(CTQMC_SEGME+CTQMC_DETSI) = op%stats(CTQMC_SEGME+CTQMC_DETSI) + SIGN(1.d0,det_ratio)
!!      END IF
!!      IF ( (time1 * beta * time_avail * DEXP(op%mu(op%Impurity%activeFlavor)*length+overlap)) &
!!           .LT. (tail * det_ratio ) ) THEN
!!        write(*,*) position
!!        CALL ImpurityOperator_remove(op%Impurity,position)
!!        CALL BathOperator_setMRemove(op%Bath,op%Impurity%particles(op%Impurity%activeFlavor))
!!        !op%seg_removed = op%seg_removed  + 1.d0
!!        op%stats(CTQMC_SEGME+CTQMC_REMOV) = op%stats(CTQMC_SEGME+CTQMC_REMOV)  + 1.d0
!!        updated = .TRUE.
!!      END IF
!!    END IF
!!  END IF
!!END SUBROUTINE Ctqmc_trySegment
!!***
!!
!!SUBROUTINE Ctqmc_tryAntiSeg(op, updated)
!!!  TYPE(BathOperator)    , INTENT(INOUT) :: Bath 
!!!  TYPE(ImpurityOperator), INTENT(INOUT) :: Impurity 
!!  LOGICAL               , INTENT(INOUT) :: updated
!!  INTEGER                               :: position
!!  DOUBLE PRECISION                      :: action
!!  DOUBLE PRECISION                      :: beta
!!  DOUBLE PRECISION                      :: time1
!!  DOUBLE PRECISION                      :: time2
!!  DOUBLE PRECISION                      :: time_avail
!!  DOUBLE PRECISION                      :: det_ratio
!!  DOUBLE PRECISION                      :: overlap
!!  DOUBLE PRECISION                      :: length
!!  DOUBLE PRECISION                      :: tail
!!  DOUBLE PRECISION, DIMENSION(1:2)      :: CdagC_1
!!
!!  IF ( .NOT. op%set ) &
!!    CALL ERROR("Ctqmc_trySegment : QMC not set                       ")
!!
!!  beta     =  op%beta
!!  tail     =  DBLE(op%Impurity%particles(op%Impurity%activeFlavor)%tail)
!!
!!  CALL RANDOM_NUMBER(action)
!!
!!  updated = .FALSE.
!!
!!  IF ( action .LT. .5d0 ) THEN ! Ajout d'un antiseg
!!    CALL RANDOM_NUMBER(time1)
!!    time1 = time1 * beta
!!    time_avail = ImpurityOperator_getAvailableTime(op%Impurity,time1,position)
!!    IF ( time_avail .LT. 0.d0 ) THEN
!!      CALL RANDOM_NUMBER(time2)
!!      time2     = time1 - time2 * time_avail
!!!      CALL CdagC_init(CdagC_1,time2,time1)
!!      CdagC_1(Cdag_) = time2
!!      CdagC_1(C_   ) = time1
!!!      length    = CdagC_length(CdagC_1) ! /!\ length is negative
!!      length = time1 - time2
!!      det_ratio = BathOperator_getDetAdd(op%Bath,CdagC_1,position,op%Impurity%particles(op%Impurity%activeFlavor))
!!      overlap   = ImpurityOperator_getNewOverlap(op%Impurity,CdagC_1) ! OK
!!      CALL RANDOM_NUMBER(time1)
!!      IF ( det_ratio .LT. 0.d0 ) THEN
!!        det_ratio    = -det_ratio
!!!        op%anti_sign = op%anti_sign + 1.d0
!!        op%stats(CTQMC_ANTIS+CTQMC_DETSI) = op%stats(CTQMC_ANTIS+CTQMC_DETSI) + SIGN(1.d0,det_ratio)
!!      END IF
!!      IF ( (time1 * (tail + 1.d0 )) & 
!!           .LT. (beta * ABS(time_avail) * det_ratio * DEXP(op%mu(op%Impurity%activeFlavor)*length + overlap) ) ) THEN
!!        CALL ImpurityOperator_add(op%Impurity,CdagC_1,position) 
!!        !write(*,*) position
!!        CALL BathOperator_setMAdd(op%bath,op%Impurity%particles(op%Impurity%activeFlavor)) 
!!        !op%anti_added = op%anti_added  + 1.d0
!!        op%stats(CTQMC_ANTIS+CTQMC_ADDED) = op%stats(CTQMC_ANTIS+CTQMC_ADDED)  + 1.d0
!!        updated = .TRUE.
!!      END IF 
!!    END IF 
!!  ELSE ! Supprimer un antiseg
!!    IF ( tail .GT. 0.d0 ) THEN
!!      CALL RANDOM_NUMBER(time1)
!!      position = -(INT(time1 * tail) + 1)
!!      time_avail = ImpurityOperator_getAvailedTime(op%Impurity,position)!OK
!!      det_ratio  = BathOperator_getDetRemove(op%Bath,position)!OK
!!      CdagC_1    = ImpurityOperator_getSegment(op%Impurity,position)!OK 
!!!      length     = CdagC_length(CdagC_1) ! /!\ length is negative
!!      length = CdagC_1(C_) - CdagC_1(Cdag_)
!!      overlap    = ImpurityOperator_getNewOverlap(op%Impurity,CdagC_1) !OK
!!      CALL RANDOM_NUMBER(time1)
!!      IF ( det_ratio .LT. 0.d0 ) THEN
!!        det_ratio   = -det_ratio
!!!        op%anti_sign = op%anti_sign + 1.d0
!!        op%stats(CTQMC_ANTIS+CTQMC_DETSI) = op%stats(CTQMC_ANTIS+CTQMC_DETSI) + SIGN(1.d0,det_ratio)
!!      END IF
!!      IF ( (time1 * beta * time_avail * DEXP(op%mu(op%Impurity%activeFlavor)*length+overlap)) &
!!           .LT. (tail * det_ratio ) ) THEN
!!        CALL ImpurityOperator_remove(op%Impurity,position)
!!        !write(*,*) position
!!        CALL BathOperator_setMRemove(op%Bath,op%Impurity%particles(op%Impurity%activeFlavor))
!!        !op%anti_removed = op%anti_removed  + 1.d0
!!        op%stats(CTQMC_ANTIS+CTQMC_REMOV) = op%stats(CTQMC_ANTIS+CTQMC_REMOV)  + 1.d0
!!        updated = .TRUE.
!!      END IF
!!    END IF
!!  END IF
!!END SUBROUTINE Ctqmc_tryAntiSeg
!!***

!!****f* ABINIT/m_Ctqmc/Ctqmc_trySwap
!! NAME
!!  Ctqmc_trySwap
!!
!! FUNCTION
!!  try a global move (swap to flavors)
!!
!! COPYRIGHT
!!  Copyright (C) 2013-2014 ABINIT group (J. Bieder)
!!  This file is distributed under the terms of the
!!  GNU General Public License, see ~abinit/COPYING
!!  or http://www.gnu.org/copyleft/gpl.txt .
!!
!! INPUTS
!!  op=ctqmc
!!
!! OUTPUT
!!  flav_i=first flavor swaped
!!  flav_j=second flavor swaped
!!
!! SIDE EFFECTS
!!
!! NOTES
!!
!! PARENTS
!!  Will be filled automatically by the parent script
!!
!! CHILDREN
!!  Will be filled automatically by the parent script
!!
!! SOURCE

SUBROUTINE Ctqmc_trySwap(op,flav_i,flav_j)

!Arguments ------------------------------------

!This section has been created automatically by the script Abilint (TD).
!Do not modify the following lines by hand.
#undef ABI_FUNC
#define ABI_FUNC 'Ctqmc_trySwap'
!End of the abilint section

  TYPE(Ctqmc)           , INTENT(INOUT) :: op
!  TYPE(BathOperator)    , INTENT(INOUT) :: Bath 
!  TYPE(ImpurityOperator), INTENT(INOUT) :: Impurity 
  INTEGER               , INTENT(  OUT) :: flav_i
  INTEGER               , INTENT(  OUT) :: flav_j
!Local variables ------------------------------
  INTEGER :: flavor_i
  INTEGER :: flavor_j
  DOUBLE PRECISION :: rnd
  DOUBLE PRECISION :: lengthi
  DOUBLE PRECISION :: lengthj
  DOUBLE PRECISION :: overlapic1
  DOUBLE PRECISION :: overlapjc1
  DOUBLE PRECISION :: overlapic2
  DOUBLE PRECISION :: overlapjc2
  DOUBLE PRECISION :: detic1
  DOUBLE PRECISION :: detjc1
  DOUBLE PRECISION :: detic2
  DOUBLE PRECISION :: detjc2
  DOUBLE PRECISION :: det_ratio
  DOUBLE PRECISION :: local_ratio

  !CALL RANDOM_NUMBER(rnd)
  CALL OurRng(op%seed,rnd)
  flavor_i = NINT(rnd*DBLE(op%flavors-1.d0))+1
  !CALL RANDOM_NUMBER(rnd)
  CALL OurRng(op%seed,rnd)
  flavor_j = NINT(rnd*DBLE(op%flavors-1.d0))+1
  
  flav_i = 0
  flav_j = 0

  IF ( flavor_i .NE. flavor_j ) THEN
    ! On tente d'intervertir i et j
    ! Configuration actuelle :
    op%modGlobalMove(2) = op%modGlobalMove(2)+1
    detic1     = BathOperator_getDetF(op%Bath,flavor_i)
    detjc1     = BathOperator_getDetF(op%Bath,flavor_j)
    lengthi    = ImpurityOperator_measN(op%Impurity,flavor_i)
    lengthj    = ImpurityOperator_measN(op%Impurity,flavor_j)
    overlapic1 = ImpurityOperator_overlapFlavor(op%Impurity,flavor_i)
    overlapjc1 = ImpurityOperator_overlapFlavor(op%Impurity,flavor_j)
    ! Configuration nouvelle :
    detic2     = BathOperator_getDetF(op%Bath,flavor_i,op%Impurity%particles(flavor_j))
    detjc2     = BathOperator_getDetF(op%Bath,flavor_j,op%Impurity%particles(flavor_i))
    ! lengths unchanged
    overlapic2 = ImpurityOperator_overlapSwap(op%Impurity,flavor_i,flavor_j)
    overlapjc2 = ImpurityOperator_overlapSwap(op%Impurity,flavor_j,flavor_i)

!    IF ( detic1*detjc1 .EQ. detic2*detjc2 ) THEN
!      det_ratio = 1.d0
!    ELSE IF ( detic1*detjc1 .EQ. 0.d0 ) THEN
!      det_ratio = detic2*detjc2 ! evite de diviser par 0 si pas de segment
!    ELSE
      det_ratio = detic2*detjc2/(detic1*detjc1)
!    END IF
    local_ratio = DEXP(-overlapic2*overlapjc2+overlapic1*overlapjc1 &
                      +(lengthj-lengthi)*(op%mu(flavor_i)-op%mu(flavor_j)))
    ! Wloc = exp(muN-Uo)
    !CALL RANDOM_NUMBER(rnd)
    CALL OurRng(op%seed,rnd)
    IF ( rnd .LT. local_ratio*det_ratio ) THEN ! swap accepted
      CALL ImpurityOperator_swap(op%Impurity, flavor_i,flavor_j)
      CALL BathOperator_swap    (op%Bath    , flavor_i,flavor_j)
      op%swap = op%swap + 1.d0
      flav_i = flavor_i
      flav_j = flavor_j
!    ELSE
!      CALL WARN("Swap refused")
!      WRITE(op%ostream,'(6E24.14)') local_ratio, det_ratio, detic1, detjc1, detic2, detjc2
    END IF
  END IF

END SUBROUTINE Ctqmc_trySwap
!!***

!!****f* ABINIT/m_Ctqmc/Ctqmc_measN
!! NAME
!!  Ctqmc_measN
!!
!! FUNCTION
!!  measure the number of electron
!!
!! COPYRIGHT
!!  Copyright (C) 2013-2014 ABINIT group (J. Bieder)
!!  This file is distributed under the terms of the
!!  GNU General Public License, see ~abinit/COPYING
!!  or http://www.gnu.org/copyleft/gpl.txt .
!!
!! INPUTS
!!  op=ctqmc
!!  iflavor=which flavor to measure
!!  updated=something has changed since last time
!!
!! OUTPUT
!!
!! SIDE EFFECTS
!!
!! NOTES
!!
!! PARENTS
!!  Will be filled automatically by the parent script
!!
!! CHILDREN
!!  Will be filled automatically by the parent script
!!
!! SOURCE

SUBROUTINE Ctqmc_measN(op, iflavor, updated)

!Arguments ------------------------------------

!This section has been created automatically by the script Abilint (TD).
!Do not modify the following lines by hand.
#undef ABI_FUNC
#define ABI_FUNC 'Ctqmc_measN'
!End of the abilint section

  TYPE(Ctqmc)             , INTENT(INOUT)     :: op
  !TYPE(ImpurityOperator), INTENT(IN   )     :: impurity
  INTEGER               , INTENT(IN   )     :: iflavor
  LOGICAL               , INTENT(IN   )     :: updated

!  IF ( .NOT. op%set ) &
!    CALL ERROR("Ctqmc_measN : QMC not set                           ")

  
  IF ( updated .EQV. .TRUE. ) THEN
    op%measN(1,iflavor) = op%measN(1,iflavor) + op%measN(3,iflavor)*op%measN(4,iflavor)
    op%measN(2,iflavor) = op%measN(2,iflavor) + op%measN(4,iflavor)
    op%measN(3,iflavor) = ImpurityOperator_measN(op%impurity)
    op%measN(4,iflavor) = 1.d0
  ELSE
    op%measN(4,iflavor) = op%measN(4,iflavor) + 1.d0
  END IF
END SUBROUTINE Ctqmc_measN
!!***

!#ifdef CTCtqmc_ANALYSIS
!!****f* ABINIT/m_Ctqmc/Ctqmc_measCorrelation
!! NAME
!!  Ctqmc_measCorrelation
!!
!! FUNCTION
!!  measure all correlations in times for a flavor
!!
!! COPYRIGHT
!!  Copyright (C) 2013-2014 ABINIT group (J. Bieder)
!!  This file is distributed under the terms of the
!!  GNU General Public License, see ~abinit/COPYING
!!  or http://www.gnu.org/copyleft/gpl.txt .
!!
!! INPUTS
!!  op=ctqmc
!!  iflavor=the flavor to measure
!!
!! OUTPUT
!!
!! SIDE EFFECTS
!!
!! NOTES
!!
!! PARENTS
!!  Will be filled automatically by the parent script
!!
!! CHILDREN
!!  Will be filled automatically by the parent script
!!
!! SOURCE

SUBROUTINE Ctqmc_measCorrelation(op, iflavor)

!Arguments ------------------------------------

!This section has been created automatically by the script Abilint (TD).
!Do not modify the following lines by hand.
#undef ABI_FUNC
#define ABI_FUNC 'Ctqmc_measCorrelation'
!End of the abilint section

  TYPE(Ctqmc)             , INTENT(INOUT)       :: op
  !TYPE(ImpurityOperator), INTENT(IN   )       :: impurity
  INTEGER               , INTENT(IN   )       :: iflavor
!Local variables ------------------------------
  INTEGER                                     :: iCdag
  INTEGER                                     :: iCdagBeta
  INTEGER                                     :: iC
  INTEGER                                     :: index
  INTEGER                                     :: size
  DOUBLE PRECISION                            :: tC
  DOUBLE PRECISION                            :: tCdag
  !DOUBLE PRECISION                            :: time
  DOUBLE PRECISION                            :: inv_dt
  DOUBLE PRECISION                            :: beta

  IF ( .NOT. op%set ) &
    CALL ERROR("Ctqmc_measCorrelation : QMC not set                 ")

  size = op%impurity%particles(op%impurity%activeFlavor)%tail
  beta = op%beta

  IF ( size .EQ. 0 ) RETURN
  
  inv_dt = op%inv_dt

  DO iCdag = 1, size ! first segments
    tCdag  = op%impurity%particles(op%impurity%activeFlavor)%list(iCdag,Cdag_)
    tC     = op%impurity%particles(op%impurity%activeFlavor)%list(iCdag,C_   )
    index = INT( ( (tC - tCdag)  * inv_dt ) + .5d0 ) + 1
    op%measCorrelation(index,1,iflavor) = op%measCorrelation(index,1,iflavor) + 1.d0
    MODCYCLE(iCdag+1,size,iCdagBeta)
    index = INT( ( ( &
                    op%impurity%particles(op%impurity%activeFlavor)%list(iCdagBeta,Cdag_) - tC &
                    + AINT(DBLE(iCdag)/DBLE(size))*beta &
                   )  * inv_dt ) + .5d0 ) + 1
    IF ( index .LT. 1 .OR. index .GT. op%samples+1 ) THEN
      CALL WARN("Ctqmc_measCorrelation : bad index line 1095         ")
    ELSE
      op%measCorrelation(index,2,iflavor) = op%measCorrelation(index,2,iflavor) + 1.d0
    END IF
!    DO iC = 1, size
!      tC = impurity%particles(impurity%activeFlavor)%list(C_,iC)
!      time = tC - tCdag
!      IF ( time .LT. 0.d0 ) time = time + beta
!      index = INT( ( time * inv_dt ) + .5d0 ) + 1
!      op%measCorrelation(index,3,iflavor) = op%measCorrelation(index,3,iflavor) + 1.d0
!    END DO
    DO iC = 1, size!  op%Greens(iflavor)%index_old%tail 
        op%measCorrelation(op%Greens(iflavor)%map%listINT(iC+(iCdag-1)*size),3,iflavor) = &
        op%measCorrelation(op%Greens(iflavor)%map%listINT(iC+(iCdag-1)*size),3,iflavor) + 1.d0
    END DO
  END DO

END SUBROUTINE Ctqmc_measCorrelation
!!***

!!****f* ABINIT/m_Ctqmc/Ctqmc_measPerturbation
!! NAME
!!  Ctqmc_measPerturbation
!!
!! FUNCTION
!!  measure perturbation order
!!
!! COPYRIGHT
!!  Copyright (C) 2013-2014 ABINIT group (J. Bieder)
!!  This file is distributed under the terms of the
!!  GNU General Public License, see ~abinit/COPYING
!!  or http://www.gnu.org/copyleft/gpl.txt .
!!
!! INPUTS
!!  op=ctqmc
!!  iflavor=the flavor to measure
!!
!! OUTPUT
!!
!! SIDE EFFECTS
!!
!! NOTES
!!
!! PARENTS
!!  Will be filled automatically by the parent script
!!
!! CHILDREN
!!  Will be filled automatically by the parent script
!!
!! SOURCE

SUBROUTINE Ctqmc_measPerturbation(op, iflavor)

!Arguments ------------------------------------

!This section has been created automatically by the script Abilint (TD).
!Do not modify the following lines by hand.
#undef ABI_FUNC
#define ABI_FUNC 'Ctqmc_measPerturbation'
!End of the abilint section

  TYPE(Ctqmc)             , INTENT(INOUT)     :: op
  !TYPE(ImpurityOperator), INTENT(IN   )     :: impurity
  INTEGER               , INTENT(IN   )     :: iflavor
!Local variables ------------------------------
  INTEGER                                   :: index

  IF ( .NOT. op%set ) &
    CALL ERROR("Ctqmc_measiPerturbation : QMC not set               ")

  index = op%impurity%particles(op%impurity%activeFlavor)%tail + 1
  IF ( index .LE. op%opt_order ) &
    op%measPerturbation(index,iflavor) = op%measPerturbation(index,iflavor) + 1.d0

END SUBROUTINE Ctqmc_measPerturbation
!!***
!#endif

!!****f* ABINIT/m_Ctqmc/Ctqmc_getResult
!! NAME
!!  Ctqmc_getResult
!!
!! FUNCTION
!!  reduce everything to get the result of the simulation
!!
!! COPYRIGHT
!!  Copyright (C) 2013-2014 ABINIT group (J. Bieder)
!!  This file is distributed under the terms of the
!!  GNU General Public License, see ~abinit/COPYING
!!  or http://www.gnu.org/copyleft/gpl.txt .
!!
!! INPUTS
!!  op=ctqmc
!!
!! OUTPUT
!!
!! SIDE EFFECTS
!!
!! NOTES
!!
!! PARENTS
!!  Will be filled automatically by the parent script
!!
!! CHILDREN
!!  Will be filled automatically by the parent script
!!
!! SOURCE

SUBROUTINE Ctqmc_getResult(op)


!This section has been created automatically by the script Abilint (TD).
!Do not modify the following lines by hand.
#undef ABI_FUNC
#define ABI_FUNC 'Ctqmc_getResult'
!End of the abilint section


#ifdef HAVE_MPI1
include 'mpif.h'
#endif
!Arguments ------------------------------------
  TYPE(Ctqmc)  , INTENT(INOUT)                    :: op
!Local variables ------------------------------
  INTEGER                                       :: iflavor
  INTEGER                                       :: flavors
  INTEGER                                       :: itau
  INTEGER                                       :: endDensity
  DOUBLE PRECISION                              :: inv_flavors
  DOUBLE PRECISION                              :: a
  DOUBLE PRECISION                              :: b
  DOUBLE PRECISION                              :: r
  DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:) :: alpha
  DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:) :: beta
  DOUBLE PRECISION,              DIMENSION(1:2) :: TabX
  DOUBLE PRECISION,              DIMENSION(1:2) :: TabY
  DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:)   :: freqs
  INTEGER, ALLOCATABLE, DIMENSION(:)   :: counts
  INTEGER, ALLOCATABLE, DIMENSION(:)   :: displs
  INTEGER                                       :: sp1
  INTEGER                                       :: spAll
  INTEGER                                       :: last
  INTEGER                                       :: n1
  INTEGER                                       :: n2
  INTEGER                                       :: debut
!  INTEGER                                       :: fin
#ifdef HAVE_MPI
  INTEGER                                       :: ierr
#endif
  DOUBLE PRECISION                              :: inv_size
  DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:) :: buffer 
  TYPE(FFTHyb) :: FFTmrka

  IF ( .NOT. op%done ) &
    CALL ERROR("Ctqmc_getResult : Simulation not run                ")

  flavors     =  op%flavors
  inv_flavors = 1.d0 / DBLE(flavors)


  inv_size = 1.d0 / DBLE(op%size)
  sp1 = 0
  spAll = 0

!#ifdef CTCtqmc_CHECK
  IF ( op%opt_check .GT. 0 ) THEN
    op%errorImpurity = ImpurityOperator_getError(op%Impurity) * inv_flavors 
    op%errorBath     = BathOperator_getError    (op%Bath    ) * inv_flavors 
  END IF
!#endif

  MALLOC(alpha,(1,1))
  MALLOC(beta,(1,1))
  MALLOC(buffer,(1,1))
  IF ( op%opt_noise .EQ. 1) THEN
    FREEIF(alpha)
    MALLOC(alpha,(1:op%samples+1,1:flavors))
    FREEIF(beta)
    MALLOC(beta,(1:op%samples+1,1:flavors))
  END IF

  IF ( op%have_MPI .EQV. .TRUE.) THEN 
    sp1   = op%samples+1
    spALL = sp1 + flavors + 6 

!#ifdef CTCtqmc_ANALYSIS
    IF ( op%opt_analysis .EQ. 1 ) &
      spAll = spAll + 3*sp1 
    IF ( op%opt_order .GT. 0 ) &
      spAll = spAll + op%opt_order 
    IF ( op%opt_noise .EQ. 1 ) &
      spAll = spAll + 2*(op%samples + 1)
!#endif

    FREEIF(buffer)
    MALLOC(buffer,(1:spAll,1:MAX(2,flavors)))
  END IF

!  op%seg_added    = op%seg_added    * inv_flavors 
!  op%seg_removed  = op%seg_removed  * inv_flavors
!  op%seg_sign     = op%seg_sign     * inv_flavors
!  op%anti_added   = op%anti_added   * inv_flavors
!  op%anti_removed = op%anti_removed * inv_flavors
!  op%anti_sign    = op%anti_sign    * inv_flavors
  op%stats(:) = op%stats(:) * inv_flavors

  DO iflavor = 1, flavors
    CALL GreenHyb_measHybrid(op%Greens(iflavor), op%Bath%M(iflavor), op%Impurity%Particles(iflavor), .TRUE.)
    CALL GreenHyb_getHybrid(op%Greens(iflavor))
    ! Accumule les dernieres mesure de N
    op%measN(1,iflavor) = op%measN(1,iflavor) + op%measN(3,iflavor)*op%measN(4,iflavor)
    op%measN(2,iflavor) = op%measN(2,iflavor) + op%measN(4,iflavor)
    ! Reduction
    op%measN(1,iflavor)  = op%measN(1,iflavor) / ( op%measN(2,iflavor) * op%beta )
    ! Correction
    CALL GreenHyb_setN(op%Greens(iflavor), op%measN(1,iflavor))
!#ifdef CTCtqmc_ANALYSIS
    IF ( op%opt_order .GT. 0 ) &
      op%measPerturbation(:   ,iflavor) = op%measPerturbation(:,iflavor) &
                                    / SUM(op%measPerturbation(:,iflavor))
    IF ( op%opt_analysis .EQ. 1 ) THEN
      op%measCorrelation (:,1,iflavor) = op%measCorrelation  (:,1,iflavor) &
                                    / SUM(op%measCorrelation (:,1,iflavor)) &
                                    * op%inv_dt 
      op%measCorrelation (:,2,iflavor) = op%measCorrelation  (:,2,iflavor) &
                                    / SUM(op%measCorrelation (:,2,iflavor)) &
                                    * op%inv_dt 
      op%measCorrelation (:,3,iflavor) = op%measCorrelation  (:,3,iflavor) &
                                    / SUM(op%measCorrelation (:,3,iflavor)) &
                                    * op%inv_dt 
    END IF
!#endif
    IF ( op%opt_noise .EQ. 1 ) THEN
      TabX(1) = DBLE(op%modNoise2)
      TabX(2) = DBLE(op%modNoise1)
      DO itau = 1, op%samples+1
        op%measNoiseG(itau,iflavor,2)%vec = -op%measNoiseG(itau,iflavor,2)%vec*op%inv_dt &  
                                           /(op%beta*DBLE(op%modNoise2))
        op%measNoiseG(itau,iflavor,1)%vec = -op%measNoiseG(itau,iflavor,1)%vec*op%inv_dt &  
                                           /(op%beta*DBLE(op%modNoise1))
        n2 = op%measNoiseG(itau,iflavor,2)%tail
        TabY(1) = Stat_deviation(op%measNoiseG(itau,iflavor,2)%vec(1:n2))!*SQRT(n2/(n2-1))
        n1 = op%measNoiseG(itau,iflavor,1)%tail
        TabY(2) = Stat_deviation(op%measNoiseG(itau,iflavor,1)%vec(1:n1))!*SQRT(n1/(n1-1))
        CALL Stat_powerReg(TabX,SQRT(2.d0*LOG(2.d0))*TabY,alpha(itau,iflavor),beta(itau,iflavor),r)
        ! ecart type -> 60%
        ! largeur a mi-hauteur d'une gaussienne -> sqrt(2*ln(2))*sigma
      END DO
    END IF

    IF ( op%have_MPI .EQV. .TRUE. ) THEN 
      buffer(1:sp1, iflavor) = op%Greens(iflavor)%oper(1:sp1)
    END IF
  END DO
  last = sp1

  op%measDE(:,:) = op%measDE(:,:) * DBLE(op%measurements) /(DBLE(op%sweeps)*op%beta)

  n1 = op%measNoise(1)%tail
  n2 = op%measNoise(2)%tail

  ! On utilise freqs comme tableau de regroupement
  ! Gather de Noise1
  IF ( op%have_MPI .EQV. .TRUE. ) THEN
    MALLOC(counts,(1:op%size))
    MALLOC(displs,(1:op%size))
    FREEIF(freqs)
    MALLOC(freqs,(1:op%size*n1))
    freqs = 0.d0
    freqs(n1*op%rank+1:n1*(op%rank+1)) = op%measNoise(1)%vec(1:n1) 
    counts(:) = n1
    displs(:) = (/ ( iflavor*n1, iflavor=0, op%size-1 ) /)
#ifdef HAVE_MPI
    CALL MPI_ALLGATHERV(MPI_IN_PLACE, 0, MPI_DOUBLE_PRECISION, &
                        freqs, counts, displs, &
                        MPI_DOUBLE_PRECISION, op%MY_COMM, ierr)
#endif
    n1 = op%size*n1
    CALL Vector_setSize(op%measNoise(1),n1)
    op%measNoise(1)%vec(1:n1) = freqs(:)
    ! Gather de Noise2
    FREE(freqs)
    MALLOC(freqs,(1:op%size*n2))
    freqs = 0.d0
    freqs(n2*op%rank+1:n2*(op%rank+1)) = op%measNoise(2)%vec(1:n2) 
    counts(:) = n2
    displs(:) = (/ ( iflavor*n2, iflavor=0, op%size-1 ) /)
#ifdef HAVE_MPI
    CALL MPI_ALLGATHERV(MPI_IN_PLACE, 0, MPI_DATATYPE_NULL, &
                        freqs, counts, displs, &
                        MPI_DOUBLE_PRECISION, op%MY_COMM, ierr)
#endif
    n2 = op%size*n2
    CALL Vector_setSize(op%measNoise(2),n2)
    op%measNoise(2)%vec(1:n2) = freqs(:)
    FREE(counts)
    FREE(displs)
    FREE(freqs)
  END IF
  !n1 = op%measNoise(1)%tail
  !n2 = op%measNoise(2)%tail

  ! Transformation des paquets pour que ca fit a CTQMC_SLICE(1|2)
  IF ( n1 .GT. CTQMC_SLICE1 ) THEN
    itau = n1/CTQMC_SLICE1
    MALLOC(freqs,(1:n1/itau))
    DO debut=1, n1/itau
      freqs(debut)=SUM(op%measNoise(1)%vec((debut-1)*itau+1:itau*debut))
    END DO
    freqs(:) = freqs(:)/DBLE(itau)
    op%modNoise1 = op%modNoise1*itau
    n1 = n1/itau
    CALL Vector_setSize(op%measNoise(1),n1)
    op%measNoise(1)%vec(1:n1) = freqs(:)
    FREE(freqs)
  END IF
  IF ( n2 .GT. CTQMC_SLICE1*CTQMC_SLICE2 ) THEN
    itau = n2/(CTQMC_SLICE1*CTQMC_SLICE2)
    MALLOC(freqs,(1:n2/itau))
    DO debut=1, n2/itau
      freqs(debut)=SUM(op%measNoise(2)%vec((debut-1)*itau+1:itau*debut))
    END DO
    freqs(:) = freqs(:)/DBLE(itau)
    op%modNoise2 = op%modNoise2*itau
    n2 = n2/itau
    CALL Vector_setSize(op%measNoise(2),n2)
    op%measNoise(2)%vec(1:n2) = freqs(:)
    FREE(freqs)
  END IF
  ! On peut s'amuser avec nos valeur d'energies
  !MALLOC(TabX,(1:20))
  !MALLOC(TabY,(1:20))

  TabX(1) = DBLE(op%modNoise2)
  TabX(2) = DBLE(op%modNoise1)

  ! Il faut calculer pour chaque modulo 10 ecarts type sur les donnes acquises
  op%measNoise(1)%vec(1:n1) = op%measNoise(1)%vec(1:n1)/(op%beta*DBLE(op%modNoise1))*DBLE(op%measurements)
  op%measNoise(2)%vec(1:n2) = op%measNoise(2)%vec(1:n2)/(op%beta*DBLE(op%modNoise2))*DBLE(op%measurements)
!  CALL Vector_print(op%measNoise(1),op%rank+70)
!  CALL Vector_print(op%measNoise(2),op%rank+50)
!  DO iflavor=1,10
!    debut = (iflavor-1)*n2/10+1
!    fin   = iflavor*n2/10
!    TabY(iflavor) = Stat_deviation(op%measNoise(2)%vec(debut:fin))
!    debut = (iflavor-1)*n1/10+1
!    fin   = iflavor*n1/10
!    TabY(10+iflavor) = Stat_deviation(op%measNoise(1)%vec(debut:fin))
!  END DO
!!  TabY(1:n) = (op%measNoise(2)%vec(1:n)   &
!!              )
!!             !/(op%beta*DBLE(op%modNoise2))*DBLE(op%measurements) &
!!             !- op%measDE(1,1))
!!  TabY(op%measNoise(2)%tail+1:n+op%measNoise(2)%tail) = (op%measNoise(1)%vec(1:n)   &
!!               )
!!             ! /(op%beta*DBLE(op%modNoise1))*DBLE(op%measurements) &
!!             ! - op%measDE(1,1))
!  IF ( op%rank .EQ. 0 ) THEN
!    DO iflavor=1,20
!      write(45,*) TabX(iflavor), TabY(iflavor)
!    END DO
!  END IF
!


  TabY(1) = Stat_deviation(op%measNoise(2)%vec(1:n2))!*SQRT(n2/(n2-1))
!!  write(op%rank+10,*) TabX(2)
!!  write(op%rank+40,*) TabX(1)
!!  CALL Vector_print(op%measNoise(1),op%rank+10)
!!  CALL Vector_print(op%measNoise(2),op%rank+40)
!!  CLOSE(op%rank+10)
!!  CLOSE(op%rank+40)
  TabY(2) = Stat_deviation(op%measNoise(1)%vec(1:n1))!*SQRT(n1/(n1-1))
!!  ! Ecart carre moyen ~ ecart type mais non biaise. Serait moins precis. Aucun
  ! impact sur la pente, juste sur l'ordonnee a l'origine.

  CALL Stat_powerReg(TabX,SQRT(2.d0*LOG(2.d0))*TabY,a,b,r)
!  FREE(TabX)
!  FREE(TabY)
  ! ecart type -> 60%
  ! largeur a mi-hauteur d'une gaussienne -> sqrt(2*ln(2))*sigma

  !op%measDE(1,1) = SUM(op%measNoise(1)%vec(1:op%measNoise(1)%tail))/(DBLE(op%measNoise(1)%tail*op%modNoise1)*op%beta)
  !op%measDE(2:flavors,1:flavors) = op%measDE(2:flavors,1:flavors) /(DBLE(op%sweeps)*op%beta)
  CALL ImpurityOperator_getErrorOverlap(op%Impurity,op%measDE)
  ! Add the difference between true calculation and quick calculation of the
  ! last sweep overlap to measDE(2,2)
  !op%measDE = op%measDE * DBLE(op%measurements) 
  IF ( op%have_MPI .EQV. .TRUE. ) THEN 
    IF ( op%opt_analysis .EQ. 1 ) THEN
      buffer(last+1:last+sp1,:) = op%measCorrelation(:,1,:)
      last = last + sp1
      buffer(last+1:last+sp1,:) = op%measCorrelation(:,2,:)
      last = last + sp1
      buffer(last+1:last+sp1,:) = op%measCorrelation(:,3,:)
      last = last + sp1
    END IF
    IF ( op%opt_order .GT. 0 ) THEN
      buffer(last+1:last+op%opt_order, :) = op%measPerturbation(:,:)
      last = last + op%opt_order
    END IF
    IF ( op%opt_noise .EQ. 1 ) THEN
      buffer(last+1:last+op%samples+1,:) = alpha(:,:)
      last = last + op%samples + 1
      buffer(last+1:last+op%samples+1,:) = beta(:,:)
      last = last + op%samples + 1
    END IF
!  op%measDE(2,2) = a*EXP(b*LOG(DBLE(op%sweeps*op%size)))
    buffer(spall-(flavors+5):spAll-6,:) = op%measDE(:,:)
!    buffer(spAll  ,1) = op%seg_added   
!    buffer(spAll-1,1) = op%seg_removed 
!    buffer(spAll-2,1) = op%seg_sign    
!    buffer(spAll  ,2) = op%anti_added  
!    buffer(spAll-1,2) = op%anti_removed
!    buffer(spAll-2,2) = op%anti_sign   
    buffer(spAll  ,1) = op%stats(1)
    buffer(spAll-1,1) = op%stats(2)
    buffer(spAll-2,1) = op%stats(3)
    buffer(spAll  ,2) = op%stats(4)
    buffer(spAll-1,2) = op%stats(5)
    buffer(spAll-2,2) = op%stats(6)
    buffer(spAll-3,1) = op%swap
    buffer(spAll-3,2) = DBLE(op%modGlobalMove(2))
    buffer(spAll-4,1) = a
    buffer(spAll-4,2) = b
!#ifdef CTCtqmc_CHECK
    buffer(spAll-5,1) = op%errorImpurity
    buffer(spAll-5,2) = op%errorBath 
!#endif

#ifdef HAVE_MPI
    CALL MPI_ALLREDUCE(MPI_IN_PLACE, buffer, spAll*flavors, &
                     MPI_DOUBLE_PRECISION, MPI_SUM, op%MY_COMM, ierr)
    CALL MPI_ALLREDUCE(MPI_IN_PLACE, op%runTime, 1, MPI_DOUBLE_PRECISION, MPI_MAX, &
             op%MY_COMM, ierr)
#endif

  
    buffer          = buffer * inv_size
    op%measDE(:,:)  = buffer(spall-(flavors+5):spAll-6,:)
!    op%seg_added    = buffer(spAll  ,1)
!    op%seg_removed  = buffer(spAll-1,1)
!    op%seg_sign     = buffer(spAll-2,1)
!    op%anti_added   = buffer(spAll  ,2)
!    op%anti_removed = buffer(spAll-1,2)
!    op%anti_sign    = buffer(spAll-2,2)
    op%stats(1)    = buffer(spAll  ,1)
    op%stats(2)    = buffer(spAll-1,1)
    op%stats(3)    = buffer(spAll-2,1)
    op%stats(4)    = buffer(spAll  ,2)
    op%stats(5)    = buffer(spAll-1,2)
    op%stats(6)    = buffer(spAll-2,2)
    op%swap         = buffer(spAll-3,1)
    op%modGlobalMove(2) = NINT(buffer(spAll-3,2))
    a               = buffer(spAll-4,1) 
    b               = buffer(spAll-4,2)
!!#ifdef CTCtqmc_CHECK
    op%errorImpurity= buffer(spAll-5,1) 
    op%errorBath    = buffer(spAll-5,2)   
!#endif

    DO iflavor = 1, flavors
      op%Greens(iflavor)%oper          = buffer(1:sp1          , iflavor)
    END DO
    last = sp1
    IF ( op%opt_analysis .EQ. 1 ) THEN
      op%measCorrelation(:,1,:) = buffer(last+1:last+sp1,:) 
      last = last + sp1
      op%measCorrelation(:,2,:) = buffer(last+1:last+sp1,:) 
      last = last + sp1
      op%measCorrelation(:,3,:) = buffer(last+1:last+sp1,:) 
      last = last + sp1
    END IF
    IF ( op%opt_order .GT. 0 ) THEN
      op%measPerturbation(:,:) = buffer(last+1:last+op%opt_order, :)
      last = last + op%opt_order
    END IF
    IF ( op%opt_noise .EQ. 1 ) THEN
      alpha(:,:) = buffer(last+1:last+op%samples+1,:)
      last = last + op%samples + 1
      beta(:,:) = buffer(last+1:last+op%samples+1,:)
      last = last + op%samples + 1
    END IF
  END IF
  DO iflavor = 1, flavors
    ! complete DE matrix
    op%measDE(iflavor, iflavor+1:flavors) = op%measDE(iflavor+1:flavors,iflavor)
  END DO
  FREE(buffer)

  IF ( op%opt_spectra .GE. 1 ) THEN
    endDensity = SIZE(op%density,2)
    IF ( op%density(1,endDensity) .EQ. -1.d0 ) &
      endDensity = endDensity - 1
    CALL FFTHyb_init(FFTmrka,endDensity,DBLE(op%thermalization)/DBLE(op%measurements*op%opt_spectra))
    ! Not very Beauty 
    MALLOC(freqs,(1:FFTmrka%size/2))
    DO iflavor = 1, flavors
      ! mean value is removed to supress the continue composent 
      CALL FFTHyb_setData(FFTmrka,op%density(iflavor,1:endDensity)/op%beta+op%Greens(iflavor)%oper(op%samples+1))
      CALL FFTHyb_run(FFTmrka,1)
      CALL FFTHyb_getData(FFTmrka,endDensity,op%density(iflavor,:),freqs)
    END DO
    op%density(flavors+1,:) = -1.d0
    op%density(flavors+1,1:FFTmrka%size/2) = freqs
    CALL FFTHyb_destroy(FFTmrka)
    FREE(freqs)
  END IF

  op%a_Noise = a
  op%b_Noise = b
  IF ( op%opt_noise .EQ. 1 ) THEN
    op%abNoiseG(1,:,:) = alpha
    op%abNoiseG(2,:,:) = beta
  END IF
  FREE(alpha)
  FREE(beta)

END SUBROUTINE Ctqmc_getResult
!!***

!!****f* ABINIT/m_Ctqmc/Ctqmc_symmetrizeGreen
!! NAME
!!  Ctqmc_symmetrizeGreen
!!
!! FUNCTION
!!  optionnaly symmetrize the green functions
!!
!! COPYRIGHT
!!  Copyright (C) 2013-2014 ABINIT group (J. Bieder)
!!  This file is distributed under the terms of the
!!  GNU General Public License, see ~abinit/COPYING
!!  or http://www.gnu.org/copyleft/gpl.txt .
!!
!! INPUTS
!!  op=ctqmc
!!  syms=weight factors
!!
!! OUTPUT
!!  argout(sizeout)=description
!!
!! SIDE EFFECTS
!!
!! NOTES
!!
!! PARENTS
!!  Will be filled automatically by the parent script
!!
!! CHILDREN
!!  Will be filled automatically by the parent script
!!
!! SOURCE

SUBROUTINE Ctqmc_symmetrizeGreen(op, syms)

!Arguments ------------------------------------

!This section has been created automatically by the script Abilint (TD).
!Do not modify the following lines by hand.
#undef ABI_FUNC
#define ABI_FUNC 'Ctqmc_symmetrizeGreen'
!End of the abilint section

  TYPE(Ctqmc)                     , INTENT(INOUT) :: op
  DOUBLE PRECISION, DIMENSION(:,:), INTENT(IN   ) :: syms
!Local variables ------------------------------
  INTEGER :: iflavor1
  INTEGER :: iflavor2
  INTEGER :: flavors
  DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:) :: green_tmp
  DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:  ) :: n_tmp

  flavors = op%flavors
  IF ( SIZE(syms,1) .NE. flavors .OR. SIZE(syms,2) .NE. flavors ) THEN
    CALL WARNALL("Ctqmc_symmetrizeGreen : wrong opt_sym -> not symmetrizing")
    RETURN
  END IF
 
  MALLOC(green_tmp,(1:op%samples+1,flavors))
  green_tmp(:,:) = 0.d0
  MALLOC(n_tmp,(1:flavors))
  n_tmp(:) = 0.d0
  DO iflavor1=1, flavors
    DO iflavor2=1,flavors
      green_tmp(:,iflavor1) = green_tmp(:,iflavor1) &
                             + syms(iflavor2,iflavor1) * op%Greens(iflavor2)%oper(:)
      n_tmp(iflavor1) = n_tmp(iflavor1) &
                             + syms(iflavor2,iflavor1) * op%measN(1,iflavor2)
    END DO
  END DO
  DO iflavor1=1, flavors
    op%Greens(iflavor1)%oper(:) = green_tmp(:,iflavor1)
    op%measN(1,iflavor1)          = n_tmp(iflavor1)
  END DO
  FREE(green_tmp)
  FREE(n_tmp)
END SUBROUTINE Ctqmc_symmetrizeGreen
!!***

!!****f* ABINIT/m_Ctqmc/Ctqmc_getGreen
!! NAME
!!  Ctqmc_getGreen
!!
!! FUNCTION
!!  Get the full green functions in time and/or frequency
!!
!! COPYRIGHT
!!  Copyright (C) 2013-2014 ABINIT group (J. Bieder)
!!  This file is distributed under the terms of the
!!  GNU General Public License, see ~abinit/COPYING
!!  or http://www.gnu.org/copyleft/gpl.txt .
!!
!! INPUTS
!!  op=ctqmc
!!
!! OUTPUT
!!  Gtau=green function in time
!!  Gw=green function in frequency
!!
!! SIDE EFFECTS
!!
!! NOTES
!!
!! PARENTS
!!  Will be filled automatically by the parent script
!!
!! CHILDREN
!!  Will be filled automatically by the parent script
!!
!! SOURCE

SUBROUTINE Ctqmc_getGreen(op, Gtau, Gw)

!Arguments ------------------------------------

!This section has been created automatically by the script Abilint (TD).
!Do not modify the following lines by hand.
#undef ABI_FUNC
#define ABI_FUNC 'Ctqmc_getGreen'
!End of the abilint section

  TYPE(Ctqmc)          , INTENT(INOUT)    :: op
  DOUBLE PRECISION, DIMENSION(:,:), OPTIONAL, INTENT(INOUT) :: Gtau
  COMPLEX(KIND=8), DIMENSION(:,:), OPTIONAL, INTENT(INOUT) :: Gw
!Local variables ------------------------------
  !INTEGER                            :: itime
  INTEGER                            :: iflavor1
  INTEGER                            :: iflavor2
  INTEGER                            :: iflavor3
  INTEGER                            :: flavors
  DOUBLE PRECISION :: u1 
  DOUBLE PRECISION :: u2
  DOUBLE PRECISION :: Un
  DOUBLE PRECISION :: UUnn

  flavors = op%flavors
  DO iflavor1 = 1, flavors
    u1 = 0.d0
    u2 = 0.d0
    DO iflavor2 = 1, flavors
      IF ( iflavor2 .EQ. iflavor1 ) CYCLE
      Un = op%Impurity%mat_U(iflavor2,iflavor1) * op%measN(1,iflavor2)
      u1 = u1 + Un 
      u2 = u2 + Un*op%Impurity%mat_U(iflavor2,iflavor1) 
      DO iflavor3 = 1, flavors
        IF ( iflavor3 .EQ. iflavor2 .OR. iflavor3 .EQ. iflavor1 ) CYCLE
        UUnn = (op%Impurity%mat_U(iflavor2,iflavor1)*op%Impurity%mat_U(iflavor3,iflavor1)) * op%measDE(iflavor2,iflavor3) 
        u2 = u2 + UUnn 
      END DO
    END DO  

    CALL GreenHyb_setMoments(op%Greens(iflavor1),u1,u2)
    IF ( PRESENT( Gtau ) ) THEN
      Gtau(1:op%samples,iflavor1) = op%Greens(iflavor1)%oper(1:op%samples)
    END IF
       !write(6,*) "present gw", present(gw)
    IF ( PRESENT( Gw ) ) THEN
       !write(6,*) "size gw",SIZE(Gw,DIM=2) ,flavors+1 
      IF ( SIZE(Gw,DIM=2) .EQ. flavors+1 ) THEN
        CALL GreenHyb_forFourier(op%Greens(iflavor1), Gomega=Gw(:,iflavor1), omega=Gw(:,op%flavors+1))
        !IF ( op%rank .EQ. 0 ) write(20,*) Gw(:,iflavor1)
      ELSE IF ( SIZE(Gw,DIM=2) .EQ. flavors ) THEN  
        CALL GreenHyb_forFourier(op%Greens(iflavor1),Gomega=Gw(:,iflavor1))
      ELSE
        CALL WARNALL("Ctqmc_getGreen : Gw is not valid                    ")
        CALL GreenHyb_forFourier(op%Greens(iflavor1),Wmax=op%Wmax)
      END IF
    ELSE
      CALL GreenHyb_forFourier(op%Greens(iflavor1),Wmax=op%Wmax)
    END IF
  END DO
END SUBROUTINE Ctqmc_getGreen
!!***

!!****f* ABINIT/m_Ctqmc/Ctqmc_getD
!! NAME
!!  Ctqmc_getD
!!
!! FUNCTION
!!  get double occupation
!!
!! COPYRIGHT
!!  Copyright (C) 2013-2014 ABINIT group (J. Bieder)
!!  This file is distributed under the terms of the
!!  GNU General Public License, see ~abinit/COPYING
!!  or http://www.gnu.org/copyleft/gpl.txt .
!!
!! INPUTS
!!  op=ctqmc
!!
!! OUTPUT
!!  D=full double occupation
!!
!! SIDE EFFECTS
!!
!! NOTES
!!
!! PARENTS
!!  Will be filled automatically by the parent script
!!
!! CHILDREN
!!  Will be filled automatically by the parent script
!!
!! SOURCE

SUBROUTINE Ctqmc_getD(op, D)

!Arguments ------------------------------------

!This section has been created automatically by the script Abilint (TD).
!Do not modify the following lines by hand.
#undef ABI_FUNC
#define ABI_FUNC 'Ctqmc_getD'
!End of the abilint section

  TYPE(Ctqmc)       , INTENT(IN ) :: op
  DOUBLE PRECISION, DIMENSION(:,:), INTENT(OUT) :: D
!Local variables ------------------------------
  INTEGER                       :: iflavor1
  INTEGER                       :: iflavor2

  IF ( SIZE(D,1) .LT. op%flavors .OR. SIZE(D,2) .LT. op%flavors ) &
    CALL ERROR("Ctqmc_getD : Dimensions of array D are too small")

  D = 0.d0

  DO iflavor1 = 1, op%flavors
    DO iflavor2 = 1, op%flavors
      D(iflavor2,iflavor1) =  op%measDE(iflavor2,iflavor1)
    END DO
    D(iflavor1,iflavor1) = 0.d0
  END DO

END SUBROUTINE Ctqmc_getD
!!***

!!****f* ABINIT/m_Ctqmc/Ctqmc_getE
!! NAME
!!  Ctqmc_getE
!!
!! FUNCTION
!!  get interaction energy and noise on it
!!
!! COPYRIGHT
!!  Copyright (C) 2013-2014 ABINIT group (J. Bieder)
!!  This file is distributed under the terms of the
!!  GNU General Public License, see ~abinit/COPYING
!!  or http://www.gnu.org/copyleft/gpl.txt .
!!
!! INPUTS
!!  op=ctqmc
!!
!! OUTPUT
!!  E=interaction energy
!!  noise=noise on this value
!!
!! SIDE EFFECTS
!!
!! NOTES
!!
!! PARENTS
!!  Will be filled automatically by the parent script
!!
!! CHILDREN
!!  Will be filled automatically by the parent script
!!
!! SOURCE

SUBROUTINE Ctqmc_getE(op,E,noise)

!Arguments ------------------------------------

!This section has been created automatically by the script Abilint (TD).
!Do not modify the following lines by hand.
#undef ABI_FUNC
#define ABI_FUNC 'Ctqmc_getE'
!End of the abilint section

  TYPE(Ctqmc)       , INTENT(IN ) :: op
  DOUBLE PRECISION, OPTIONAL, INTENT(OUT) :: E
  DOUBLE PRECISION, OPTIONAL, INTENT(OUT) :: Noise

  IF ( PRESENT(E) ) &
    E = op%measDE(1,1)  
  IF ( PRESENT(Noise) ) &
    Noise = SUM(op%Impurity%mat_U)/(op%flavors*(op%flavors-1)) &
            * op%a_Noise*(DBLE(op%sweeps)*DBLE(op%size))**op%b_Noise
END SUBROUTINE Ctqmc_getE
!!***

!!****f* ABINIT/m_Ctqmc/Ctqmc_printAll
!! NAME
!!  Ctqmc_printAll
!!
!! FUNCTION
!!  print different functions computed during the simulation
!!
!! COPYRIGHT
!!  Copyright (C) 2013-2014 ABINIT group (J. Bieder)
!!  This file is distributed under the terms of the
!!  GNU General Public License, see ~abinit/COPYING
!!  or http://www.gnu.org/copyleft/gpl.txt .
!!
!! INPUTS
!!  op=ctqmc
!!
!! OUTPUT
!!
!! SIDE EFFECTS
!!
!! NOTES
!!
!! PARENTS
!!  Will be filled automatically by the parent script
!!
!! CHILDREN
!!  Will be filled automatically by the parent script
!!
!! SOURCE

SUBROUTINE Ctqmc_printAll(op)

!Arguments ------------------------------------

!This section has been created automatically by the script Abilint (TD).
!Do not modify the following lines by hand.
#undef ABI_FUNC
#define ABI_FUNC 'Ctqmc_printAll'
!End of the abilint section

  TYPE(Ctqmc), INTENT(INOUT) :: op

  IF ( .NOT. op%done ) &
    CALL WARNALL("Ctqmc_printAll : Simulation not run                 ")

  CALL Ctqmc_printQMC(op)

  CALL Ctqmc_printGreen(op)

  CALL Ctqmc_printD(op)

!  CALL Ctqmc_printE(op)

!#ifdef CTCtqmc_ANALYSIS
  CALL Ctqmc_printPerturbation(op)

  CALL Ctqmc_printCorrelation(op)
!#endif

  CALL Ctqmc_printSpectra(op)

END SUBROUTINE Ctqmc_printAll
!!***

!!****f* ABINIT/m_Ctqmc/Ctqmc_printQMC
!! NAME
!!  Ctqmc_printQMC
!!
!! FUNCTION
!!  print ctqmc statistics
!!
!! COPYRIGHT
!!  Copyright (C) 2013-2014 ABINIT group (J. Bieder)
!!  This file is distributed under the terms of the
!!  GNU General Public License, see ~abinit/COPYING
!!  or http://www.gnu.org/copyleft/gpl.txt .
!!
!! INPUTS
!!  op=ctqmc
!!
!! OUTPUT
!!
!! SIDE EFFECTS
!!
!! NOTES
!!
!! PARENTS
!!  Will be filled automatically by the parent script
!!
!! CHILDREN
!!  Will be filled automatically by the parent script
!!
!! SOURCE

SUBROUTINE Ctqmc_printQMC(op)

!Arguments ------------------------------------

!This section has been created automatically by the script Abilint (TD).
!Do not modify the following lines by hand.
#undef ABI_FUNC
#define ABI_FUNC 'Ctqmc_printQMC'
!End of the abilint section

  TYPE(Ctqmc), INTENT(INOUT) :: op
!Local variables ------------------------------
  INTEGER                  :: ostream
  INTEGER                  :: iflavor
  DOUBLE PRECISION         :: sweeps
  DOUBLE PRECISION         :: invSweeps
  CHARACTER(LEN=2)         :: a
  CHARACTER(LEN=15)        :: string

  !IF ( op%rank .NE. 0) RETURN
  IF ( op%rank .NE. MOD(op%size,op%size)) RETURN

  ostream   = op%ostream
  sweeps    = DBLE(op%sweeps)
  invSweeps = 1.d0/sweeps

  WRITE(ostream,'(1x,F13.0,A11,F10.2,A12,I5,A5)') sweeps*DBLE(op%size), " sweeps in ", op%runTime, &
                 " seconds on ", op%size, " CPUs"
  WRITE(ostream,'(A28,F6.2)') "Segments added        [%] : ", op%stats(4)*invSweeps*100.d0
  WRITE(ostream,'(A28,F6.2)') "Segments removed      [%] : ", op%stats(5)*invSweeps*100.d0
  WRITE(ostream,'(A28,F6.2)') "Segments sign         [%] : ", op%stats(6)*invSweeps*100.d0
  WRITE(ostream,'(A28,F6.2)') "Anti-segments added   [%] : ", op%stats(1)*invSweeps*100.d0
  WRITE(ostream,'(A28,F6.2)') "Anti-segments removed [%] : ", op%stats(2)*invSweeps*100.d0
  WRITE(ostream,'(A28,F6.2)') "Anti-segments sign    [%] : ", op%stats(3)*invSweeps*100.d0
  IF ( op%modGlobalMove(1) .LT. op%sweeps + 1 ) THEN
    WRITE(ostream,'(A28,F6.2)') "Global Move           [%] : ", op%swap         *invSweeps*100.d0*op%modGlobalMove(1)
    WRITE(ostream,'(A28,F6.2)') "Global Move Reduced   [%] : ", op%swap         / DBLE(op%modGlobalMove(2))*100.d0
  END IF
!#ifdef CTCtqmc_CHECK
  IF ( op%opt_check .EQ. 1 .OR. op%opt_check .EQ. 3 ) &
    WRITE(ostream,'(A28,E22.14)') "Impurity test         [%] : ", op%errorImpurity*100.d0
  IF ( op%opt_check .GE. 2 ) &
      WRITE(ostream,'(A28,E22.14)') "Bath     test         [%] : ", op%errorBath    *100.d0
!#endif
  WRITE(ostream,'(A28,ES22.14,A5,ES21.14)') "<Epot>                [U] : ", op%measDE(1,1), " +/- ",&
!#ifdef HAVE_MPI
         SUM(op%Impurity%mat_U)/(op%flavors*(op%flavors-1)) * op%a_Noise*(sweeps*DBLE(op%size))**op%b_Noise
!#else
!                                                              op%a_Noise*(sweeps)**op%b_Noise
!#endif
  WRITE(ostream,'(A28,F8.4,A3,F7.4)') "Noise                [/U] : ", op%a_Noise, " x^", op%b_Noise
  WRITE(ostream,'(A28,E10.2)')  "Niquist puls.     [/beta] : ", ACOS(-1.d0)*op%inv_dt
  WRITE(ostream,'(A28,E22.14)') "Max Acc. Epot Error   [U] : ", op%measDE(2,2)/(op%beta*op%modNoise1*2.d0)*sweeps
  
  !WRITE(ostream,'(A28,F7.4,A3,F7.4,A4,E20.14)') "Noise            [G(tau)] : ", op%a_Noise(2), "x^", op%b_Noise(2), " -> ", &
                                                              !op%a_Noise(2)*(sweeps*DBLE(op%size))**op%b_Noise(2)
  IF ( op%opt_order .GT. 0 ) THEN 
    WRITE(a,'(I2)') op%flavors
    string = '(A28,'//TRIM(ADJUSTL(a))//'(1x,I3))'
    WRITE(ostream,string) "Perturbation orders       : ",(/ (MAXLOC(op%measPerturbation(:, iflavor))-1, iflavor=1, op%flavors) /)
  END IF
  !CALL FLUSH(op%ostream)
  IF ( ABS(((op%stats(4) *invSweeps*100.d0) / (op%stats(5) *invSweeps*100.d0) - 1.d0)) .GE. 0.02d0 &
   .OR. ABS(((op%stats(1)*invSweeps*100.d0) / (op%stats(2)*invSweeps*100.d0) - 1.d0)) .GE. 0.02d0 ) &
    THEN 
    CALL WARNALL("Ctqmc_printQMC : bad statistic according to moves. Increase sweeps")
  END IF
  IF ( ABS(op%b_Noise+0.5)/0.5d0 .GE. 0.05d0 ) &
    CALL WARNALL("Ctqmc_printQMC : bad statistic according to Noise. Increase sweeps")
!  IF ( ISNAN(op%a_Noise) .OR. ISNAN(op%a_Noise) ) &
!    CALL WARNALL("Ctqmc_printQMC : NaN appeared. Increase sweeps    ")


END SUBROUTINE Ctqmc_printQMC
!!***

!!****f* ABINIT/m_Ctqmc/Ctqmc_printGreen
!! NAME
!!  Ctqmc_printGreen
!!
!! FUNCTION
!!  print green functions
!!
!! COPYRIGHT
!!  Copyright (C) 2013-2014 ABINIT group (J. Bieder)
!!  This file is distributed under the terms of the
!!  GNU General Public License, see ~abinit/COPYING
!!  or http://www.gnu.org/copyleft/gpl.txt .
!!
!! INPUTS
!!  op=ctqmc
!!  oFileIn=file stream
!!
!! OUTPUT
!!
!! SIDE EFFECTS
!!
!! NOTES
!!
!! PARENTS
!!  Will be filled automatically by the parent script
!!
!! CHILDREN
!!  Will be filled automatically by the parent script
!!
!! SOURCE

SUBROUTINE Ctqmc_printGreen(op, oFileIn)

!Arguments ------------------------------------

!This section has been created automatically by the script Abilint (TD).
!Do not modify the following lines by hand.
#undef ABI_FUNC
#define ABI_FUNC 'Ctqmc_printGreen'
!End of the abilint section

  TYPE(Ctqmc)        , INTENT(IN)    :: op
  INTEGER  , OPTIONAL, INTENT(IN)    :: oFileIn
!Local variables ------------------------------
  INTEGER                            :: oFile
  INTEGER                            :: itime
  INTEGER                            :: sp1
  INTEGER                            :: iflavor
  INTEGER                            :: flavors
  CHARACTER(LEN=4)                   :: cflavors
  CHARACTER(LEN=50)                  :: string
  DOUBLE PRECISION                   :: dt
  DOUBLE PRECISION                   :: sweeps

  !IF ( op%rank .NE. MOD(1,op%size)) RETURN
  IF ( op%rank .NE. MOD(op%size+1,op%size)) RETURN

  oFile = 40
  IF ( PRESENT(oFileIn) ) THEN
    oFile = oFileIn
  ELSE
    OPEN(UNIT=oFile, FILE="Gtau.dat")
  END IF

  sp1     =  op%samples
  dt      =  op%beta / DBLE(sp1)
  sp1     =  sp1 + 1
  flavors =  op%flavors
  sweeps = DBLE(op%sweeps)*DBLE(op%size)

  IF ( op%opt_noise .EQ. 1) THEN
    WRITE(cflavors,'(I4)') 2*flavors+1
    string = '(1x,'//TRIM(ADJUSTL(cflavors))//'ES22.14)'
    DO itime = 1, sp1
      WRITE(oFile,string) DBLE(itime-1)*dt, &
      (/ (op%Greens(iflavor)%oper(itime), iflavor=1, flavors) /), &
      (/ (op%abNoiseG(1,itime,iflavor)*(sweeps)**op%abNoiseG(2,itime,iflavor), iflavor=1, flavors) /)
    END DO
  ELSE
    WRITE(cflavors,'(I4)') flavors+1
    string = '(1x,'//TRIM(ADJUSTL(cflavors))//'ES22.14)'
    DO itime = 1, sp1
      WRITE(oFile,string) DBLE(itime-1)*dt, &
      (/ (op%Greens(iflavor)%oper(itime), iflavor=1, flavors) /)
    END DO
  END IF

  IF ( .NOT. PRESENT(oFileIn) ) CLOSE(oFile)

END SUBROUTINE Ctqmc_printGreen
!!***

!!****f* ABINIT/m_Ctqmc/Ctqmc_printD
!! NAME
!!  Ctqmc_printD
!!
!! FUNCTION
!!  print individual double occupancy
!!
!! COPYRIGHT
!!  Copyright (C) 2013-2014 ABINIT group (J. Bieder)
!!  This file is distributed under the terms of the
!!  GNU General Public License, see ~abinit/COPYING
!!  or http://www.gnu.org/copyleft/gpl.txt .
!!
!! INPUTS
!!  op=ctqmc
!!  oFileIn=file stream
!!
!! OUTPUT
!!
!! SIDE EFFECTS
!!
!! NOTES
!!
!! PARENTS
!!  Will be filled automatically by the parent script
!!
!! CHILDREN
!!  Will be filled automatically by the parent script
!!
!! SOURCE

SUBROUTINE Ctqmc_printD(op,oFileIn)

!Arguments ------------------------------------

!This section has been created automatically by the script Abilint (TD).
!Do not modify the following lines by hand.
#undef ABI_FUNC
#define ABI_FUNC 'Ctqmc_printD'
!End of the abilint section

  TYPE(Ctqmc)          , INTENT(IN)    :: op
  INTEGER  , OPTIONAL, INTENT(IN)    :: oFileIn
!Local variables ------------------------------
  INTEGER                            :: oFile
  INTEGER                            :: iflavor1
  INTEGER                            :: iflavor2

  !IF ( op%rank .NE. MOD(2,op%size)) RETURN
  IF ( op%rank .NE. MOD(op%size+2,op%size)) RETURN

  oFile = 41
  IF ( PRESENT(oFileIn) ) THEN
    oFile = oFileIn
  ELSE
    OPEN(UNIT=oFile, FILE="D.dat")
  END IF

  DO iflavor1 = 1, op%flavors
    DO iflavor2 = iflavor1+1, op%flavors
      WRITE(oFile,'(1x,A8,I4,A1,I4,A3,ES21.14)') "Orbitals", iflavor1, "-", iflavor2, " : ", op%measDE(iflavor2,iflavor1)
    END DO
  END DO

  IF ( .NOT. PRESENT(oFileIn) ) CLOSE(oFile)

END SUBROUTINE Ctqmc_printD
!!***

!!****f* ABINIT/m_Ctqmc/Ctqmc_printE
!! NAME
!!  Ctqmc_printE
!!
!! FUNCTION
!!  print energy and noise 
!!
!! COPYRIGHT
!!  Copyright (C) 2013-2014 ABINIT group (J. Bieder)
!!  This file is distributed under the terms of the
!!  GNU General Public License, see ~abinit/COPYING
!!  or http://www.gnu.org/copyleft/gpl.txt .
!!
!! INPUTS
!!  op=ctqmc
!!  oFileIn=file stream
!!
!! OUTPUT
!!
!! SIDE EFFECTS
!!
!! NOTES
!!
!! PARENTS
!!  Will be filled automatically by the parent script
!!
!! CHILDREN
!!  Will be filled automatically by the parent script
!!
!! SOURCE

SUBROUTINE Ctqmc_printE(op,oFileIn)

!Arguments ------------------------------------

!This section has been created automatically by the script Abilint (TD).
!Do not modify the following lines by hand.
#undef ABI_FUNC
#define ABI_FUNC 'Ctqmc_printE'
!End of the abilint section

  TYPE(Ctqmc)          , INTENT(IN)    :: op
  INTEGER  , OPTIONAL, INTENT(IN)    :: oFileIn
!Local variables ------------------------------
  INTEGER                            :: oFile
  DOUBLE PRECISION                   :: E
  DOUBLE PRECISION                   :: Noise

  !IF ( op%rank .NE. MOD(3,op%size)) RETURN
  IF ( op%rank .NE. MOD(op%size+3,op%size)) RETURN

  oFile = 42
  IF ( PRESENT(oFileIn) ) THEN
    oFile = oFileIn
  ELSE
    OPEN(UNIT=oFile, FILE="BetaENoise.dat")
  END IF

  CALL Ctqmc_getE(op,E,Noise)

  WRITE(oFile,'(1x,F5.2,A2,ES21.14,A2,ES21.14)') op%beta, "  ", E, "  ",  Noise

  IF ( .NOT. PRESENT(oFileIn) ) CLOSE(oFile)

END SUBROUTINE Ctqmc_printE
!!***

!#ifdef CTCtqmc_ANALYSIS
!!****f* ABINIT/m_Ctqmc/Ctqmc_printPerturbation
!! NAME
!!  Ctqmc_printPerturbation
!!
!! FUNCTION
!!  print perturbation order
!!
!! COPYRIGHT
!!  Copyright (C) 2013-2014 ABINIT group (J. Bieder)
!!  This file is distributed under the terms of the
!!  GNU General Public License, see ~abinit/COPYING
!!  or http://www.gnu.org/copyleft/gpl.txt .
!!
!! INPUTS
!!  op=ctqmc
!!  oFileIn=file stream
!!
!! OUTPUT
!!  argout(sizeout)=description
!!
!! SIDE EFFECTS
!!
!! NOTES
!!
!! PARENTS
!!  Will be filled automatically by the parent script
!!
!! CHILDREN
!!  Will be filled automatically by the parent script
!!
!! SOURCE

SUBROUTINE Ctqmc_printPerturbation(op, oFileIn)

!Arguments ------------------------------------

!This section has been created automatically by the script Abilint (TD).
!Do not modify the following lines by hand.
#undef ABI_FUNC
#define ABI_FUNC 'Ctqmc_printPerturbation'
!End of the abilint section

  TYPE(Ctqmc)          , INTENT(IN)           :: op
  INTEGER  , OPTIONAL,  INTENT(IN)          :: oFileIn
!Local variables-------------------------------
  INTEGER                                   :: oFile
  INTEGER                                   :: iorder
  INTEGER                                   :: order
  INTEGER                                   :: iflavor
  INTEGER                                   :: flavors
  CHARACTER(LEN=2)                          :: a
  CHARACTER(LEN=50)                         :: string

  !IF ( op%rank .NE. MOD(4,op%size)) RETURN
  IF ( op%rank .NE. MOD(op%size+4,op%size)) RETURN
  IF ( op%opt_order .LE. 0 ) RETURN

  oFile = 43
  IF ( PRESENT(oFileIn) ) THEN
    oFile = oFileIn
  ELSE
    OPEN(UNIT=oFile, FILE="Perturbation.dat")
  END IF
    
  order        =  op%opt_order
  flavors      =  op%flavors

  WRITE(a,'(I2)') flavors
  string = '(I5,'//TRIM(ADJUSTL(a))//'F19.15)'
  DO iorder = 1, order
    WRITE(oFile,string) iorder-1, &
                (/ (op%measPerturbation(iorder, iflavor), iflavor=1, flavors) /)
  END DO

  IF ( .NOT. PRESENT(oFileIn) ) CLOSE(oFile)
END SUBROUTINE Ctqmc_printPerturbation
!!***

!!****f* ABINIT/m_Ctqmc/Ctqmc_printCorrelation
!! NAME
!!  Ctqmc_printCorrelation
!!
!! FUNCTION
!!  print correlation fonctions
!!
!! COPYRIGHT
!!  Copyright (C) 2013-2014 ABINIT group (J. Bieder)
!!  This file is distributed under the terms of the
!!  GNU General Public License, see ~abinit/COPYING
!!  or http://www.gnu.org/copyleft/gpl.txt .
!!
!! INPUTS
!!  op=ctqmc
!!  oFileIn=file stream
!!
!! OUTPUT
!!
!! SIDE EFFECTS
!!
!! NOTES
!!
!! PARENTS
!!  Will be filled automatically by the parent script
!!
!! CHILDREN
!!  Will be filled automatically by the parent script
!!
!! SOURCE

SUBROUTINE Ctqmc_printCorrelation(op, oFileIn)

!Arguments ------------------------------------

!This section has been created automatically by the script Abilint (TD).
!Do not modify the following lines by hand.
#undef ABI_FUNC
#define ABI_FUNC 'Ctqmc_printCorrelation'
!End of the abilint section

  TYPE(Ctqmc)          , INTENT(IN)             :: op
  INTEGER  , OPTIONAL, INTENT(IN)             :: oFileIn
!Local variables ------------------------------
  INTEGER                                     :: oFile
  INTEGER                                     :: itime
  INTEGER                                     :: sp1
  INTEGER                                     :: iflavor
  INTEGER                                     :: i
  INTEGER                                     :: flavors
  CHARACTER(LEN=2)                            :: a
  CHARACTER(LEN=50)                           :: string
  DOUBLE PRECISION                            :: dt

  !IF ( op%rank .NE. MOD(5,op%size)) RETURN
  IF ( op%rank .NE. MOD(op%size+5,op%size)) RETURN
  IF ( op%opt_analysis .NE. 1 ) RETURN

  oFile = 44
  IF ( PRESENT(oFileIn) ) THEN
    oFile = oFileIn
  ELSE
    OPEN(UNIT=oFile, FILE="Correlation.dat")
  END IF

  sp1         =  op%samples
  dt          =  op%beta / sp1
  sp1         =  sp1 + 1
  flavors     =  op%flavors

  i = 3*flavors + 1
  WRITE(a,'(I2)') i
  WRITE(oFile,*) "# time  (/ (segement, antiseg, correl), i=1, flavor/)"
  string = '(1x,'//TRIM(ADJUSTL(a))//'F19.15)'
  DO itime = 1, sp1
    WRITE(oFile,string) DBLE(itime-1)*dt, &
                   (/ ( &
                   (/ ( op%measCorrelation(itime, i, iflavor), i=1,3) /) &
                   , iflavor=1, flavors) /)
  END DO

  IF ( .NOT. PRESENT(oFileIn) ) CLOSE(oFile)

END SUBROUTINE Ctqmc_printCorrelation
!!***
!#endif

!!****f* ABINIT/m_Ctqmc/Ctqmc_printSpectra
!! NAME
!!  Ctqmc_printSpectra
!!
!! FUNCTION
!!  print fourier transform of time evolution of number of electrons
!!
!! COPYRIGHT
!!  Copyright (C) 2013-2014 ABINIT group (J. Bieder)
!!  This file is distributed under the terms of the
!!  GNU General Public License, see ~abinit/COPYING
!!  or http://www.gnu.org/copyleft/gpl.txt .
!!
!! INPUTS
!!  op=ctqmc
!!  oFileIn=file stream
!!
!! OUTPUT
!!
!! SIDE EFFECTS
!!
!! NOTES
!!
!! PARENTS
!!  Will be filled automatically by the parent script
!!
!! CHILDREN
!!  Will be filled automatically by the parent script
!!
!! SOURCE

SUBROUTINE Ctqmc_printSpectra(op, oFileIn)

!Arguments ------------------------------------

!This section has been created automatically by the script Abilint (TD).
!Do not modify the following lines by hand.
#undef ABI_FUNC
#define ABI_FUNC 'Ctqmc_printSpectra'
!End of the abilint section

  TYPE(Ctqmc)          , INTENT(IN)             :: op
  INTEGER  , OPTIONAL, INTENT(IN)             :: oFileIn
!Local variables ------------------------------
  INTEGER                                     :: oFile
  INTEGER                                     :: flavors
  INTEGER                                     :: indDensity
  INTEGER                                     :: endDensity
  CHARACTER(LEN=4)                            :: a
  CHARACTER(LEN=16)                           :: formatSpectra

  !IF ( op%rank .NE. MOD(6,op%size)) RETURN
  IF ( op%opt_spectra .LT. 1 ) RETURN

  oFile = 45+op%rank
  a ="0000"
  WRITE(a,'(I4)') op%rank
  IF ( PRESENT(oFileIn) ) THEN
    oFile = oFileIn
  ELSE
    OPEN(UNIT=oFile, FILE="Markov_"//TRIM(ADJUSTL(a))//".dat")
  END IF

  flavors     =  op%flavors
  WRITE(a,'(I4)') flavors+1
  formatSpectra ='(1x,'//TRIM(ADJUSTL(a))//'ES22.14)'
  WRITE(oFile,*) "# freq[/hermalization] FFT"

  endDensity = SIZE(op%density,2)
  DO WHILE ( op%density(flavors+1,endDensity) .EQ. -1 )
    endDensity = endDensity -1
  END DO

  DO indDensity = 1, endDensity
    WRITE(oFile,formatSpectra) op%density(flavors+1,indDensity), op%density(1:flavors,indDensity)
  END DO

  IF ( .NOT. PRESENT(oFileIn) ) CLOSE(oFile)

END SUBROUTINE Ctqmc_printSpectra
!!***

!!****f* ABINIT/m_Ctqmc/Ctqmc_destroy
!! NAME
!!  Ctqmc_destroy
!!
!! FUNCTION
!!  destroy and deallocate all variables
!!
!! COPYRIGHT
!!  Copyright (C) 2013-2014 ABINIT group (J. Bieder)
!!  This file is distributed under the terms of the
!!  GNU General Public License, see ~abinit/COPYING
!!  or http://www.gnu.org/copyleft/gpl.txt .
!!
!! INPUTS
!!  op=ctqmc
!!
!! OUTPUT
!!
!! SIDE EFFECTS
!!
!! NOTES
!!
!! PARENTS
!!  Will be filled automatically by the parent script
!!
!! CHILDREN
!!  Will be filled automatically by the parent script
!!
!! SOURCE

SUBROUTINE Ctqmc_destroy(op)

!Arguments ------------------------------------

!This section has been created automatically by the script Abilint (TD).
!Do not modify the following lines by hand.
#undef ABI_FUNC
#define ABI_FUNC 'Ctqmc_destroy'
!End of the abilint section

  TYPE(Ctqmc), INTENT(INOUT) :: op
!Local variables ------------------------------
  INTEGER                  :: iflavor
  INTEGER                  :: flavors
  INTEGER                  :: i
  INTEGER                  :: j
  INTEGER                  :: k

  if ( op%init .EQV. .FALSE. ) RETURN

  flavors = op%flavors

  CALL ImpurityOperator_destroy(op%Impurity)
  CALL BathOperator_destroy(op%Bath)
  CALL Vector_destroy(op%measNoise(1))
  CALL Vector_destroy(op%measNoise(2))

  IF ( ALLOCATED(op%Greens) ) THEN
    DO iflavor = 1, flavors
     CALL GreenHyb_destroy(op%Greens(iflavor))
    END DO
    DT_FREE( op%Greens )
  END IF
!#ifdef CTCtqmc_ANALYSIS
  FREEIF(op%measCorrelation)
  FREEIF(op%measPerturbation)
  FREEIF(op%measN)
  FREEIF(op%measDE)
  FREEIF(op%mu)
  FREEIF(op%abNoiseG)
  IF ( ALLOCATED(op%measNoiseG) ) THEN
    DO i=1,2
      DO j = 1, op%flavors
        DO k= 1, op%samples+1
          CALL Vector_destroy(op%measNoiseG(k,j,i))
        END DO
      END DO
    END DO
    DT_FREE(op%measNoiseG)
  END IF
  FREEIF(op%density)
!#endif
  op%ostream        = 0
  op%istream        = 0
 
  op%sweeps         = 0
  op%thermalization = 0
  op%flavors        = 0
  op%samples        = 0
  op%beta           = 0.d0
!  op%seg_added      = 0.d0
!  op%anti_added     = 0.d0
!  op%seg_removed    = 0.d0
!  op%anti_removed   = 0.d0
!  op%seg_sign       = 0.d0
!  op%anti_sign      = 0.d0
  op%stats          = 0.d0
  op%swap           = 0.d0


  op%set  = .FALSE.
  op%done = .FALSE.
  op%init = .FALSE.
END SUBROUTINE Ctqmc_destroy
!!***

END MODULE m_Ctqmc
!!***
