#if defined HAVE_CONFIG_H
#include "config.h"
#endif
!{\src2tex{textfont=tt}}
!!****m* ABINIT/m_ImpurityOperator
!! NAME
!!  m_ImpurityOperator
!! 
!! FUNCTION 
!!  manage all related to Impurity
!!
!! 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_ImpurityOperator
USE m_ListCdagC
USE m_Global
IMPLICIT NONE

!!***

PRIVATE

!!****t* m_ImpurityOperator/ImpurityOperator
!! NAME
!!  ImpurityOperator
!!
!! 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 :: ImpurityOperator
  LOGICAL          _PRIVATE :: doCheck = .FALSE.
  INTEGER          _PRIVATE :: flavors
  INTEGER                   :: activeFlavor
  DOUBLE PRECISION _PRIVATE          :: beta
  DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:)          :: mat_U 
  DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:) _PRIVATE :: overlaps   ! total overlaps
  DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:  ) _PRIVATE :: updates    ! new_(anti)seg 
  TYPE(ListCdagC)                               _PRIVATE :: list_swap
  TYPE(ListCdagC) , ALLOCATABLE, DIMENSION(:  )          :: particles 
  DOUBLE PRECISION _PRIVATE :: checkNumber
  DOUBLE PRECISION _PRIVATE :: tolerance
  DOUBLE PRECISION _PRIVATE :: meanError
END TYPE ImpurityOperator
!!***

PUBLIC  :: ImpurityOperator_init
PUBLIC  :: ImpurityOperator_reset
PUBLIC  :: ImpurityOperator_computeU
PUBLIC  :: ImpurityOperator_setUmat
PUBLIC  :: ImpurityOperator_activateParticle
PUBLIC  :: ImpurityOperator_getAvailableTime
PUBLIC  :: ImpurityOperator_getAvailedTime
PUBLIC  :: ImpurityOperator_add
PUBLIC  :: ImpurityOperator_getSegment
PUBLIC  :: ImpurityOperator_remove
PUBLIC  :: ImpurityOperator_getNewOverlap
PRIVATE :: ImpurityOperator_overlapSegFlav
PUBLIC  :: ImpurityOperator_overlapflavor
PUBLIC  :: ImpurityOperator_overlapSwap
PUBLIC  :: ImpurityOperator_swap
PRIVATE :: ImpurityOperator_overlapIJ
PUBLIC  :: ImpurityOperator_measDE
PUBLIC  :: ImpurityOperator_cleanOverlaps
PUBLIC  :: ImpurityOperator_measN
PUBLIC  :: ImpurityOperator_destroy
PUBLIC  :: ImpurityOperator_getErrorOverlap
PUBLIC  :: ImpurityOperator_doCheck
PRIVATE :: ImpurityOperator_checkOverlap
PUBLIC  :: ImpurityOperator_getError
PUBLIC  :: ImpurityOperator_printLatex

CONTAINS
!!***

!SUBROUTINE ImpurityOperator_init(op, flavors, beta, N)
!!****f* ABINIT/m_ImpurityOperator/ImpurityOperator_init
!! NAME
!!  ImpurityOperator_init
!!
!! FUNCTION
!!  Initialize and allocate
!!
!! 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=ImpurtiyOperator
!!  flavors=number of flavors
!!  beta=inverse temperature
!!
!! OUTPUT
!!
!! SIDE EFFECTS
!!
!! NOTES
!!
!! PARENTS
!!  Will be filled automatically by the parent script
!!
!! CHILDREN
!!  Will be filled automatically by the parent script
!!
!! SOURCE

SUBROUTINE ImpurityOperator_init(op, flavors, beta)

!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 'ImpurityOperator_init'
!End of the abilint section

  TYPE(ImpurityOperator), INTENT(INOUT) :: op
  INTEGER               , INTENT(IN   ) :: flavors
  !DOUBLE PRECISION      , INTENT(IN   ) :: U
  !DOUBLE PRECISION      , INTENT(IN   ) :: J
  DOUBLE PRECISION      , INTENT(IN   ) :: beta
  !INTEGER               , INTENT(IN   ) :: N
!Local variables ------------------------------
  INTEGER                               :: IT
  
  op%flavors      = flavors
  op%activeFlavor = 0
  op%beta         = beta

  IF ( MOD(flavors,2) .NE. 0 ) &
    CALL ERROR("ImpurityOperator_init : flavors is not even        ")

!#ifdef CTQMC_CHECK
  op%meanError    = 0.d0
  op%checkNumber  = 0.d0
  op%tolerance    = 0.d0
  op%doCheck      = .FALSE.
!#endif
  DT_FREEIF(op%particles)
  DT_MALLOC(op%particles,(1:flavors))
  FREEIF(op%mat_U)
  MALLOC(op%mat_U,(1:flavors,1:flavors))
  FREEIF(op%overlaps)
  MALLOC(op%overlaps,(1:flavors,1:flavors))
  op%overlaps = 0.d0
  FREEIF(op%updates)
  MALLOC(op%updates,(1:flavors))
  op%updates = 0.d0
  !CALL ImpurityOperator_computeU(op, U, J)
  !op%mat_U = U
  !IF ( ASSOCIATED(op%mu) ) FREE(op%mu)
  !MALLOC(op%mu,(1:flavors))
 
  !op%shift_mu = SUM(op%mat_U(:,1)) * .5d0 
  DO IT = 1,flavors
    !CALL ListCdagC_init(op%particles(IT), DBLE(N)/beta,100) !FIXME size of the List
    CALL ListCdagC_init(op%particles(IT),100) !FIXME size of the List
    op%particles(IT)%list(0,C_   ) = beta ! Empty orbital 
    op%particles(IT)%list(0,Cdag_) = 0.d0
!    op%particles(IT)%list(0)%Cdag = beta ! Full orbital 
!    op%particles(IT)%list(0)%C    = 0.d0
  END DO
  op%activeFlavor = 0
END SUBROUTINE ImpurityOperator_init 
!!***

!!****f* ABINIT/m_ImpurityOperator/ImpurityOperator_reset
!! NAME
!!  ImpurityOperator_reset
!!
!! FUNCTION
!!  reset operator
!!
!! 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=ImpurtiyOperator
!!
!! OUTPUT
!!
!! SIDE EFFECTS
!!
!! NOTES
!!
!! PARENTS
!!  Will be filled automatically by the parent script
!!
!! CHILDREN
!!  Will be filled automatically by the parent script
!!
!! SOURCE

SUBROUTINE ImpurityOperator_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 'ImpurityOperator_reset'
!End of the abilint section

  TYPE(ImpurityOperator), INTENT(INOUT) :: op
!Local variables ------------------------------
  INTEGER                               :: IT

  op%activeFlavor = 0
  op%overlaps = 0.d0
  op%updates = 0.d0
!#ifdef CTQMC_CHECK
  op%meanError    = 0.d0
  op%checkNumber  = 0.d0
  op%tolerance    = 0.d0
  op%doCheck      = .FALSE.
!#endif
  DO IT = 1,op%flavors
    CALL ListCdagC_clear(op%particles(IT)) 
    op%particles(IT)%list(0,C_   )    = op%beta ! Empty orbital 
    op%particles(IT)%list(0,Cdag_) = 0.d0
!    op%particles(IT)%list(0)%Cdag = beta ! Full orbital 
!    op%particles(IT)%list(0)%C    = 0.d0
  END DO

END SUBROUTINE ImpurityOperator_reset
!!***

!!****f* ABINIT/m_ImpurityOperator/ImpurityOperator_computeU
!! NAME
!!  ImpurityOperator_computeU
!!
!! FUNCTION
!!  Compute an interaction matrix for t2g like interaction
!!
!! 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=ImpurityOperator
!!  U=Coulomb scrren interaction
!!  J=Hund couplage
!!
!! OUTPUT
!!
!! SIDE EFFECTS
!!
!! NOTES
!!
!! PARENTS
!!  Will be filled automatically by the parent script
!!
!! CHILDREN
!!  Will be filled automatically by the parent script
!!
!! SOURCE

SUBROUTINE ImpurityOperator_computeU(op, U, 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 'ImpurityOperator_computeU'
!End of the abilint section

  TYPE(ImpurityOperator), INTENT(INOUT) :: op
  DOUBLE PRECISION      , INTENT(IN   ) :: U
  DOUBLE PRECISION      , INTENT(IN   ) :: J
!Local variables ------------------------------
  INTEGER                               :: flavor11
  INTEGER                               :: flavor12
  INTEGER                               :: flavor21
  INTEGER                               :: flavor22
  INTEGER                               :: flavors
  INTEGER                               :: flavors_2
  DOUBLE PRECISION                      :: Uprime

  Uprime = U - 2.d0 * J
  flavors = op%flavors
  flavors_2 = flavors / 2
  DO flavor11 = 1, flavors_2
    flavor12 = flavors - flavor11 + 1
    op%mat_U(flavor11, flavor11) = 0.d0
    op%mat_U(flavor12, flavor12) = 0.d0
    op%mat_U(flavor11+flavors_2, flavor11) = U
    op%mat_U(flavor12-flavors_2, flavor12) = U
    DO flavor21 = flavor11+1, flavors_2
      flavor22 = flavors - flavor21 + 1
      op%mat_U(flavor21, flavor11) = Uprime
      op%mat_U(flavor22-flavors_2, flavor12-flavors_2) = Uprime
      op%mat_U(flavor21+flavors_2, flavor11+flavors_2) = Uprime
      op%mat_U(flavor22, flavor12) = Uprime
    END DO
    DO flavor21 = flavor11+flavors_2+1, flavors
      flavor22 = flavors - flavor21 + 1
      op%mat_U(flavor21, flavor11) = Uprime - J
      op%mat_U(flavor22+flavors_2, flavor12-flavors_2) = Uprime - J
      op%mat_U(flavor21-flavors_2, flavor11+flavors_2) = Uprime - J
      op%mat_U(flavor22, flavor12) = Uprime - J
    END DO
  END DO
END SUBROUTINE ImpurityOperator_computeU 
!!***

!!****f* ABINIT/m_ImpurityOperator/ImpurityOperator_setUmat
!! NAME
!!  ImpurityOperator_setUmat
!!
!! FUNCTION
!!  Set directly the U 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=ImpurtityOperator
!!  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 ImpurityOperator_setUmat(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 'ImpurityOperator_setUmat'
!End of the abilint section

  TYPE(ImpurityOperator), INTENT(INOUT) :: op
  DOUBLE PRECISION, DIMENSION(:,:), INTENT(IN   ) :: matU

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

  op%mat_U = matU
END SUBROUTINE ImpurityOperator_setUmat
!!***

!!****f* ABINIT/m_ImpurityOperator/ImpurityOperator_activateParticle
!! NAME
!!  ImpurityOperator_activateParticle
!!
!! FUNCTION
!!  active 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=ImpurityOperator
!!  flavor=the flavor
!!
!! OUTPUT
!!
!! SIDE EFFECTS
!!
!! NOTES
!!
!! PARENTS
!!  Will be filled automatically by the parent script
!!
!! CHILDREN
!!  Will be filled automatically by the parent script
!!
!! SOURCE

SUBROUTINE ImpurityOperator_activateParticle(op,flavor)

!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 'ImpurityOperator_activateParticle'
!End of the abilint section

  TYPE(ImpurityOperator), INTENT(INOUT) :: op
  INTEGER               , INTENT(IN   ) :: flavor

  IF ( flavor .GT. op%flavors ) &
    CALL ERROR("ImpurityOperator_activateParticle : out of range  ")
  IF ( ALLOCATED(op%particles) ) THEN 
    op%activeFlavor   =  flavor
  ELSE
    CALL ERROR("ImpurityOperator_activateParticle : not allocated  ")
  END IF
END SUBROUTINE ImpurityOperator_activateParticle
!!***

!!****f* ABINIT/m_ImpurityOperator/ImpurityOperator_getAvailableTime
!! NAME
!!  ImpurityOperator_getAvailableTime
!!
!! FUNCTION
!!  get the time available and the position of the segment to consider
!!  negative if on a segment
!!  positive if outside a 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=ImpurityOperator
!!  time=time to look for
!!
!! OUTPUT
!!  ImpurityOperator_getAvailableTime=Time available
!!  position=position of the next segment
!!
!! SIDE EFFECTS
!!
!! NOTES
!!
!! PARENTS
!!  Will be filled automatically by the parent script
!!
!! CHILDREN
!!  Will be filled automatically by the parent script
!!
!! SOURCE

DOUBLE PRECISION FUNCTION ImpurityOperator_getAvailableTime(op, time, position)

!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 'ImpurityOperator_getAvailableTime'
!End of the abilint section

  TYPE(ImpurityOperator), INTENT(IN   ) :: op
  DOUBLE PRECISION      , INTENT(IN   ) :: time
  INTEGER               , INTENT(OUT  ) :: position
!Local variables ------------------------------
  DOUBLE PRECISION                      :: t_avail
  INTEGER                               :: position_dwn
  INTEGER                               :: aF
#include "ListCdagC_firstHigher.h"
  aF = op%activeFlavor
  IF ( aF .LE. 0 ) &
    CALL ERROR("ImpurityOperator_getAvailableTime : no active flav")
  
  IF ( op%particles(aF)%tail .EQ. 0 ) THEN
    t_avail = op%particles(aF)%list(0,C_) - op%particles(aF)%list(0,Cdag_)
    position = SIGN(1,INT(t_avail))
  ELSE
!    position = ListCdagC_firstHigher( op%particles(aF), time ) 
#define list_1 op%particles(aF) 
#include "ListCdagC_firstHigher"
#undef list_1
    position = firstHigher
    position_dwn = position - 1
    IF ( position_dwn .LE. 0) position_dwn = op%particles(aF)%tail
  
!    t_avail = (time - op%particles(aF)%list(position_dwn)) .MOD. op%beta
    t_avail = time - op%particles(aF)%list(position_dwn,C_)
    IF ( op%particles(aF)%list(position_dwn,Cdag_) .GT. time ) &
      t_avail = t_avail + op%beta 
  
    IF ( t_avail .GT. 0.d0 ) THEN  !! We are outside the position_dwn segment
!      t_avail = (op%particles(aF)%list(ABS(position)) - time ) .MOD. op%beta 
      t_avail = op%particles(aF)%list(ABS(position),Cdag_) - time 
      IF ( op%particles(aF)%list(ABS(position),Cdag_) .LT. time ) &
        t_avail = t_avail + op%beta
      ! ABS is used to prevent position to be -1 which is HERE the same as 1
    ELSE
      position = - position_dwn
    END IF
  END IF
  
    ImpurityOperator_getAvailableTime = t_avail

END FUNCTION ImpurityOperator_getAvailableTime
!!***

!!****f* ABINIT/m_ImpurityOperator/ImpurityOperator_getAvailedTime
!! NAME
!!  ImpurityOperator_getAvailedTime
!!
!! FUNCTION
!!  get the time available without the segment "position"
!!
!! 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=ImpurityOperator
!!  position=position of the segment
!!
!! OUTPUT
!!  ImpurityOperator_getAvailedTime=time available before ...
!!
!! SIDE EFFECTS
!!
!! NOTES
!!
!! PARENTS
!!  Will be filled automatically by the parent script
!!
!! CHILDREN
!!  Will be filled automatically by the parent script
!!
!! SOURCE

DOUBLE PRECISION FUNCTION ImpurityOperator_getAvailedTime(op, position)

!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 'ImpurityOperator_getAvailedTime'
!End of the abilint section

  TYPE(ImpurityOperator), INTENT(IN   ) :: op
  INTEGER               , INTENT(IN   ) :: position
  DOUBLE PRECISION                      :: T_avail
  INTEGER                               :: Pup
  INTEGER                               :: ABSp
  INTEGER                               :: tail
  INTEGER                               :: aF

  aF = op%activeFlavor
  IF ( aF .LE. 0 ) &
    CALL ERROR("ImpurityOperator_getAvailedTime : no active flavor")
  ABSp = ABS(position)
!  position_up = (ABSposition+1).MOD.op%particles(aF)%tail
 tail = op%particles(aF)%tail
  MODCYCLE(ABSp+1,tail,Pup)
  IF ( position .GT. 0 ) THEN
    t_avail = op%particles(aF)%list(Pup, Cdag_) &
            - op%particles(aF)%list(ABSp,Cdag_)
  ELSE
    t_avail = op%particles(aF)%list(Pup ,C_) &
            - op%particles(aF)%list(ABSp,C_)
  END IF
  IF ( t_avail .LE. 0.d0 ) t_avail = t_avail + op%beta
  ImpurityOperator_getAvailedTime = t_avail
END FUNCTION ImpurityOperator_getAvailedTime
!!***

!!****f* ABINIT/m_ImpurityOperator/ImpurityOperator_add
!! NAME
!!  ImpurityOperator_add
!!
!! FUNCTION
!!  add a segment to the active 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=ImpurityOperator
!!  CdagC_1=couple of times
!!  position_val=position of the CdagC_1 couple in the list
!!
!! OUTPUT
!!
!! SIDE EFFECTS
!!
!! NOTES
!!
!! PARENTS
!!  Will be filled automatically by the parent script
!!
!! CHILDREN
!!  Will be filled automatically by the parent script
!!
!! SOURCE

SUBROUTINE ImpurityOperator_add(op, CdagC_1, position_val)

!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 'ImpurityOperator_add'
!End of the abilint section

  TYPE(ImpurityOperator), INTENT(INOUT) :: op
  DOUBLE PRECISION, DIMENSION(1:2), INTENT(IN   ) :: CdagC_1
  INTEGER               , INTENT(IN   ) :: position_val
!Local variables ------------------------------
  INTEGER                               :: position
  INTEGER                               :: aF
  INTEGER                               :: i
  DOUBLE PRECISION, DIMENSION(1:2)      :: C2modify
  DOUBLE PRECISION, DIMENSION(1:2)      :: C2add
  DOUBLE PRECISION                      :: TCdag
  DOUBLE PRECISION                      :: TC

  aF = op%activeFlavor
  IF ( aF .LE. 0 ) &
    CALL ERROR("ImpurityOperator_add : no active flavor           ")
  
  position = position_val

  IF ( CdagC_1(C_) .GT. CdagC_1(Cdag_) ) THEN ! Ajout d'un segment
    C2add = CdagC_1
  ELSE                                        ! Ajout d'un antisegment
    IF ( (op%particles(aF)%tail .EQ. 0) .AND. (op%particles(aF)%list(0,C_) .EQ. 0d0)) THEN ! should be full orbital
      IF ( CdagC_1(Cdag_) .GT. op%beta ) THEN
!        CALL CdagC_init(C2add,CdagC_1%Cdag-op%beta,CdagC_1%C)
        C2add(Cdag_) = CdagC_1(Cdag_)-op%beta
        C2add(C_   ) = CdagC_1(C_)
      ELSE
!        CALL CdagC_init(C2add,CdagC_1%Cdag,CdagC_1%C+op%beta)
        C2add(Cdag_) = CdagC_1(Cdag_)
        C2add(C_   ) = CdagC_1(C_)+op%beta
      END IF
      position = 0
      op%particles(aF)%list(0,C_   ) = op%beta
      op%particles(aF)%list(0,Cdag_) = 0.d0
    ELSE IF ( op%particles(aF)%tail .GT. 0 ) THEN
      position = ABS(position)
      TCdag = op%particles(aF)%list(position,Cdag_)
      TC    = CdagC_1(C_)
      IF ( TCdag .GT. TC ) TC = TC + op%beta
!      CALL CdagC_init(C2modify,TCdag,TC)
      C2modify(Cdag_) = TCdag
      C2modify(C_   ) = TC
  
!      TCdag    = CdagC_1%Cdag.MOD.op%beta
      MODCYCLE(CdagC_1(Cdag_),op%beta,TCdag)
      TC       = op%particles(aF)%list(position,C_)
!      CALL CdagC_init(C2add,TCdag,TC)
      C2add(Cdag_) = TCdag
      C2add(C_   ) = TC
  
      op%particles(aF)%list(position,:) = C2modify
      IF ( C2modify(Cdag_) .GT. C2add(Cdag_) ) THEN
        position = 0
!        C2add%C = C2add%C.MOD.op%beta
        MODCYCLE(C2add(C_),op%beta,C2add(C_))
      END IF
    ELSE
      CALL ERROR("ImpurityOperator_add : try to add an antisegment to an empty orbital")
    END IF
    position = position + 1
  END IF
  CALL ListCdagC_insert(op%particles(aF), c2add, position)
  DO i = 1, op%flavors
    op%overlaps(i,aF) = op%overlaps(i,aF) + op%updates(i)
    op%overlaps(aF,i) = op%overlaps(i,aF)
  END DO

END SUBROUTINE ImpurityOperator_add
!!***

!!****f* ABINIT/m_ImpurityOperator/ImpurityOperator_getSegment
!! NAME
!!  ImpurityOperator_getSegment
!!
!! FUNCTION
!!  Return the segment at position_val
!!
!! 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=ImpurityOperator
!!  position_val=position of the asked segment
!!
!! OUTPUT
!!  ImpurityOperator_getSegment(2)=the couple of time
!!
!! SIDE EFFECTS
!!
!! NOTES
!!
!! PARENTS
!!  Will be filled automatically by the parent script
!!
!! CHILDREN
!!  Will be filled automatically by the parent script
!!
!! SOURCE

FUNCTION ImpurityOperator_getSegment(op,position_val)

!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 'ImpurityOperator_getSegment'
!End of the abilint section

  TYPE(ImpurityOperator), INTENT(INOUT) :: op
  INTEGER               , INTENT(IN   ) :: position_val
!Local variables ------------------------------
  INTEGER                               :: position
  INTEGER                               :: tail
  INTEGER                               :: aF
  DOUBLE PRECISION                      :: beta
  DOUBLE PRECISION                      :: ImpurityOperator_getSegment(1:2)

  aF = op%activeFlavor
  IF ( aF .LE. 0 ) &
    CALL ERROR("ImpurityOperator_getSegment : no active flavor    ")

  IF ( position_val .GT. 0 ) THEN
    ImpurityOperator_getSegment = op%particles(aF)%list(position_val,1:2)
  ELSE
    position = ABS(position_val)
    tail = op%particles(aF)%tail
    beta = op%beta
    ImpurityOperator_getSegment(C_)  = op%particles(aF)%list(position,C_)
    position = position + 1
    IF ( position .GT. tail ) THEN
      IF ( ImpurityOperator_getSegment(C_) .LT. beta ) THEN
        ImpurityOperator_getSegment(Cdag_) = op%particles(aF)%list(1,Cdag_) + beta
      ELSE
        ImpurityOperator_getSegment(Cdag_) = op%particles(aF)%list(1,Cdag_)
        ImpurityOperator_getSegment(C_)    = ImpurityOperator_getSegment(C_) -beta
      END IF
    ELSE
      ImpurityOperator_getSegment(Cdag_) = op%particles(aF)%list(position,Cdag_)
    END IF

  END IF
END FUNCTION ImpurityOperator_getSegment
!!***

!!****f* ABINIT/m_ImpurityOperator/ImpurityOperator_remove
!! NAME
!!  ImpurityOperator_remove
!!
!! FUNCTION
!!  Remove a segment for the active 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=ImpurityOperator
!!  ieme=segment to remove
!!
!! OUTPUT
!!
!! SIDE EFFECTS
!!
!! NOTES
!!
!! PARENTS
!!  Will be filled automatically by the parent script
!!
!! CHILDREN
!!  Will be filled automatically by the parent script
!!
!! SOURCE

SUBROUTINE ImpurityOperator_remove(op,ieme)

!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 'ImpurityOperator_remove'
!End of the abilint section

  TYPE(ImpurityOperator), INTENT(INOUT) :: op
  INTEGER               , INTENT(IN   ) :: ieme
!Local variables ------------------------------
  DOUBLE PRECISION, DIMENSION(1:2)      :: CdagC_1
  INTEGER                               :: position
  INTEGER                               :: position_dwn
  INTEGER                               :: i
  INTEGER                               :: tail
  INTEGER                               :: aF
!  DOUBLE PRECISION                      :: toRemove

  aF = op%activeFlavor
  IF ( aF .LE. 0 ) &
    CALL ERROR("ImpurityOperator_removeIeme : no active flavor    ")
  position = ABS(ieme)
  IF ( position .GT. op%particles(aF)%tail ) &
    CALL ERROR("ImpurityOperator_removeIeme : out of range        ")

  IF ( (ieme .LT. 0)  .AND. (op%particles(aF)%tail .GT. 1) ) THEN 
    position_dwn = position
!    position = (position+1).MOD.op%particles(aF)%tail
    tail = op%particles(aF)%tail
    MODCYCLE((position+1),tail,position)
    CdagC_1(Cdag_) = op%particles(aF)%list(position_dwn,Cdag_)
    CdagC_1(C_   ) = op%particles(aF)%list(position,C_)
    IF (position_dwn .GT. position) CdagC_1(C_) = CdagC_1(C_) + op%beta
!    toRemove  = op%particles(aF)%list(position)%C - (CdagC_1%C.MOD.op%beta)
!    CdagC_1%C = CdagC_1%C + toRemove  
    op%particles(aF)%list(position_dwn,:) = CdagC_1
  END IF

  IF ( position .EQ. 1 ) THEN
    SELECT CASE (ieme)
      CASE (1) 
        op%particles(aF)%list(0,C_   ) = op%beta
        op%particles(aF)%list(0,Cdag_) = 0.d0 
      CASE (-1)
        op%particles(aF)%list(0,C_   ) = 0.d0 
        op%particles(aF)%list(0,Cdag_) = op%beta
    END SELECT
  END IF
  CALL ListCdagC_erase(op%particles(aF),position)
  DO i = 1, op%flavors
    op%overlaps(i,aF) = op%overlaps(i,aF) - op%updates(i)
    op%overlaps(aF,i) = op%overlaps(i,aF)
  END DO
END SUBROUTINE ImpurityOperator_remove
!!***

!!****f* ABINIT/m_ImpurityOperator/ImpurityOperator_getNewOverlap
!! NAME
!!  ImpurityOperator_getNewOverlap
!!
!! FUNCTION
!!  Get the overlap induced by CdagC_1 in the current configuration
!!
!! 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=ImpurityOperator
!!  CdagC_1=the segment
!!
!! OUTPUT
!!  ImpurityOperator_getNewOverlap=overlap..
!!
!! SIDE EFFECTS
!!
!! NOTES
!!
!! PARENTS
!!  Will be filled automatically by the parent script
!!
!! CHILDREN
!!  Will be filled automatically by the parent script
!!
!! SOURCE

DOUBLE PRECISION FUNCTION ImpurityOperator_getNewOverlap(op, CdagC_1)

!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 'ImpurityOperator_getNewOverlap'
!End of the abilint section

  TYPE(ImpurityOperator), INTENT(INOUT) :: op
  DOUBLE PRECISION, DIMENSION(1:2), INTENT(IN) :: CdagC_1
!Local variables ------------------------------
  DOUBLE PRECISION, DIMENSION(1:2)   :: CdagC_2
  DOUBLE PRECISION                   :: overlap
  DOUBLE PRECISION                   :: totalOverlap
  DOUBLE PRECISION                   :: sign
  INTEGER                            :: flavor
  INTEGER                            :: otherFlavor

  flavor = op%activeFlavor
  IF ( flavor .LE. 0 ) &
    CALL ERROR("ImpurityOperator_getNewOverlap : no active flavor ")
  IF ( CdagC_1(Cdag_) .LT. CdagC_1(C_) ) THEN ! segment C*C
    CdagC_2 = CdagC_1
    sign = -1.d0
  ELSE
    CdagC_2(C_) = CdagC_1(Cdag_)
    CdagC_2(Cdag_) = CdagC_1(C_)
    sign = 1.d0
  END IF

  totalOverlap = 0.d0

  DO otherFlavor = 1, op%flavors
    IF ( otherFlavor .EQ. flavor ) CYCLE
    overlap = ImpurityOperator_overlapSegFlav(op,CdagC_2(1:2),otherflavor)
    totalOverlap = totalOverlap &
                 + overlap * op%mat_U(otherFlavor,flavor)
    op%updates(otherFlavor) = -sign * overlap
  END DO

  totalOverlap = totalOverlap * sign
  ImpurityOperator_getNewOverlap = totalOverlap

END FUNCTION ImpurityOperator_getNewOverlap
!!***

!!****f* ABINIT/m_ImpurityOperator/ImpurityOperator_overlapSegFlav
!! NAME
!!  ImpurityOperator_overlapSegFlav
!!
!! FUNCTION
!!  Compute the overlap of a segment with 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=ImpurityOperator
!!  CdagC_1=segment
!!  flavor=flavor to use
!!
!! OUTPUT
!!  ImpurityOperator_overlapSegFlav=overlap between CdagC_1 and flavor
!!
!! SIDE EFFECTS
!!
!! NOTES
!!
!! PARENTS
!!  Will be filled automatically by the parent script
!!
!! CHILDREN
!!  Will be filled automatically by the parent script
!!
!! SOURCE

DOUBLE PRECISION FUNCTION ImpurityOperator_overlapSegFlav(op,CdagC_1,flavor)

!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 'ImpurityOperator_overlapSegFlav'
!End of the abilint section

  TYPE(ImpurityOperator), INTENT(INOUT) :: op
!  TYPE(CdagC)           , INTENT(IN) :: CdagC_1
  DOUBLE PRECISION, DIMENSION(1:2), INTENT(IN) :: CdagC_1
  INTEGER               , INTENT(IN) :: flavor
!Local variables ------------------------------
!  TYPE(CdagC), DIMENSION(:), POINTER :: list => NULL()
  DOUBLE PRECISION                   :: totalCdag
  DOUBLE PRECISION                   :: totalC
  DOUBLE PRECISION                   :: beta
  DOUBLE PRECISION                   :: Time
  DOUBLE PRECISION                   :: Tmin
  DOUBLE PRECISION                   :: Tmax
  DOUBLE PRECISION                   :: TmaxBeta
  DOUBLE PRECISION                   :: TscanMin
  DOUBLE PRECISION                   :: TscanMax
!  DOUBLE PRECISION                   :: sign
  DOUBLE PRECISION                   :: C
  DOUBLE PRECISION                   :: Cdag
  DOUBLE PRECISION                   :: loop
  DOUBLE PRECISION                   :: itmin
  DOUBLE PRECISION                   :: itmax
  INTEGER                            :: tail
  INTEGER                            :: tp1
  INTEGER                            :: scanning 
  INTEGER                            :: imin 
  INTEGER                            :: imax 
  INTEGER                            :: loops
  INTEGER                            :: iloop
#include "ListCdagC_firstHigher.h"

  beta = op%beta
  Tmin = CdagC_1(Cdag_)
  Tmax = CdagC_1(C_)
  itmin= 0.d0

!  TmaxBeta     = Tmax.MOD.beta
  MODCYCLE(Tmax,beta,TmaxBeta)

  tail = op%particles(flavor)%tail

  totalC = 0.d0
  totalCdag = 0.d0
  IF ( tail .NE. 0 ) THEN
    tp1  = tail + 1
    loop = 0.d0
!    imin = ListCdagC_firstHigher( op%particles(flavor), Tmin ) - 1
    Time = Tmin
#define list_1 op%particles(flavor) 
#include "ListCdagC_firstHigher"
    imin = firstHigher - 1

    SELECT CASE ( imin ) 
      CASE(0)
        scanning = tail
        loop = -1.d0
      CASE(-2)
        scanning = tail
      CASE DEFAULT
        scanning = imin
    END SELECT
!    imax = ListCdagC_firstHigher( op%particles(flavor), TmaxBeta ) !- 1 Jamais atteint
    Time = TmaxBeta
#include "ListCdagC_firstHigher"
#undef list_1
    imax = firstHigher

    TscanMin = Tmin
    TscanMax = Tmax

    ! Regarder avant 
    IF ( (imin .EQ. 0) ) THEN
      C = op%particles(flavor)%list(scanning,C_) +loop*  beta
      Cdag = op%particles(flavor)%list(scanning,Cdag_) +loop* beta
      itmax = MAX(TscanMin, Cdag)
      itmin = MIN(TscanMax, C   )

      IF ( itmin .GT. itmax ) THEN ! si egal alors overlap de 0
        totalC = totalC + itmin
        totalCdag = totalCdag + itmax
      END IF
      scanning = scanning+1
      IF ( scanning .EQ. tp1 ) THEN
        scanning = 1
      END IF
    END IF

    loops = imax - scanning
    IF ( TmaxBeta .NE. Tmax ) THEN
      loops = tail - loops
    ELSE IF ( imax .EQ. -1 ) THEN
      loops = tail - imin
    END IF

    !Comparer betement 2 segments
    DO iloop =0, loops
      C = op%particles(flavor)%list(scanning,C_)
      Cdag = op%particles(flavor)%list(scanning,Cdag_)
      itmax = MAX(TscanMin, Cdag)
      itmin = MIN(TscanMax,C)

      IF ( itmin .GT. itmax ) THEN ! si egal alors overla de 0
        totalC = totalC + itmin
        totalCdag = totalCdag + itmax
      END IF
      scanning = scanning + 1
      IF ( scanning .EQ. tp1 ) THEN
        scanning = 1
        IF ( itmin .EQ. TScanMax ) EXIT
        TscanMin = TscanMin - beta
        TscanMax = TscanMax - beta
      END IF
    END DO

    ! Regarder apres le segment
    IF ( (itmin .NE. TscanMax) ) THEN
      C = op%particles(flavor)%list(scanning,C_)
      Cdag = op%particles(flavor)%list(scanning,Cdag_) 
      itmax = MAX(TscanMin, Cdag)
      itmin = MIN(TscanMax,C)

      IF ( itmin .GT. itmax ) THEN ! si egal alors overla de 0
        totalC = totalC + itmin
        totalCdag = totalCdag + itmax
      END IF
    END IF
  ELSE IF ( op%particles(flavor)%list(0,C_) .EQ. 0.d0 ) THEN ! full orbital
      totalC    = Tmax
      totalCdag = Tmin
  END IF
!#ifdef CTQMC_CHECK
  IF ( op%doCheck .EQV. .TRUE. ) &
    CALL ImpurityOperator_checkOverlap(op, Tmin, Tmax,totalC-totalCdag,flavor)
!#endif
  ImpurityOperator_overlapSegFlav = totalC - totalCdag 

END FUNCTION ImpurityOperator_overlapSegFlav
!!***

!!****f* ABINIT/m_ImpurityOperator/ImpurityOperator_overlapFlavor
!! NAME
!!  ImpurityOperator_overlapFlavor
!!
!! FUNCTION
!!  Returns the overlap of flavor with the others
!!
!! 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=ImpurityOperator
!!  flavor=the one we want
!!
!! OUTPUT
!!  ImpurityOperator_overlapFlavor=result
!!
!! SIDE EFFECTS
!!
!! NOTES
!!
!! PARENTS
!!  Will be filled automatically by the parent script
!!
!! CHILDREN
!!  Will be filled automatically by the parent script
!!
!! SOURCE

DOUBLE PRECISION FUNCTION ImpurityOperator_overlapFlavor(op,flavor)

!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 'ImpurityOperator_overlapFlavor'
!End of the abilint section

  TYPE(ImpurityOperator), INTENT(IN) :: op
  INTEGER,      OPTIONAL, INTENT(IN) :: flavor
!Local variables ------------------------------
  INTEGER                            :: otherFlavor
  DOUBLE PRECISION                   :: overlap
  DOUBLE PRECISION                   :: totalOverlap

  totalOverlap = 0.d0
  DO otherFlavor = 1, op%flavors
    IF ( otherFlavor .EQ. flavor ) CYCLE
    overlap = op%overlaps(otherFlavor,flavor)
    totalOverlap = totalOverlap &
                 + overlap * op%mat_U(otherFlavor,flavor)
  END DO

  ImpurityOperator_overlapFlavor = totalOverlap

END FUNCTION ImpurityOperator_overlapflavor
!!***

!!****f* ABINIT/m_ImpurityOperator/ImpurityOperator_overlapSwap
!! NAME
!!  ImpurityOperator_overlapSwap
!!
!! FUNCTION
!!  compute the overlap of flavor1 with the configuration of flavor2
!!
!! 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=ImpurityOperator
!!  flavor1=interaction value
!!  flavor2=configuration
!!
!! OUTPUT
!!  ImpurityOperator_overlapSwap=new overlap
!!
!! SIDE EFFECTS
!!
!! NOTES
!!
!! PARENTS
!!  Will be filled automatically by the parent script
!!
!! CHILDREN
!!  Will be filled automatically by the parent script
!!
!! SOURCE

DOUBLE PRECISION FUNCTION ImpurityOperator_overlapSwap(op,flavor1,flavor2)

!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 'ImpurityOperator_overlapSwap'
!End of the abilint section

  TYPE(ImpurityOperator), INTENT(IN) :: op
  INTEGER               , INTENT(IN) :: flavor1
  INTEGER               , INTENT(IN) :: flavor2
!Local variables ------------------------------
  INTEGER                            :: otherFlavor
  DOUBLE PRECISION                   :: overlap
  DOUBLE PRECISION                   :: totalOverlap

  totalOverlap = 0.d0
! Calcul l'overlap de flavor1 en utilisant la configuration de flavor2
  DO otherFlavor = 1, op%flavors
    IF ( otherFlavor .EQ. flavor2 ) THEN
      CYCLE
    ELSE IF ( otherFlavor .EQ. flavor1 ) THEN
      overlap = op%overlaps(otherFlavor,flavor2)
      totalOverlap = totalOverlap &
                   + overlap * op%mat_U(otherFlavor,flavor2)
    ELSE
      overlap = op%overlaps(otherFlavor,flavor2)
      totalOverlap = totalOverlap &
                   + overlap * op%mat_U(otherFlavor,flavor1)
    END IF
  END DO

  ImpurityOperator_overlapSwap = totalOverlap

END FUNCTION ImpurityOperator_overlapSwap
!!***

!!****f* ABINIT/m_ImpurityOperator/ImpurityOperator_swap
!! NAME
!!  ImpurityOperator_swap
!!
!! FUNCTION
!!  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=ImpurtiyOperator
!!  flavor1=to swap
!!  flavor2=to swap
!!
!! OUTPUT
!!
!! SIDE EFFECTS
!!
!! NOTES
!!
!! PARENTS
!!  Will be filled automatically by the parent script
!!
!! CHILDREN
!!  Will be filled automatically by the parent script
!!
!! SOURCE

SUBROUTINE ImpurityOperator_swap(op,flavor1, flavor2)

!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 'ImpurityOperator_swap'
!End of the abilint section

  TYPE(ImpurityOperator), INTENT(INOUT) :: op
  INTEGER               , INTENT(IN   ) :: flavor1
  INTEGER               , INTENT(IN   ) :: flavor2
!Local variables ------------------------------
  INTEGER                               :: iflavor
  DOUBLE PRECISION                      :: overlap_tmp

  DO iflavor = 1, op%flavors
    IF ( iflavor .NE. flavor1  .AND. iflavor .NE. flavor2) THEN
      overlap_tmp = op%overlaps(iflavor,flavor1)
      op%overlaps(iflavor,flavor1) = op%overlaps(iflavor,flavor2)
      op%overlaps(flavor1,iflavor) = op%overlaps(iflavor,flavor2)
      op%overlaps(iflavor,flavor2) = overlap_tmp
      op%overlaps(flavor2,iflavor) = overlap_tmp
    END IF
  END DO

  !CALL ListCdagC_print(op%particles(flavor1),233)
  !CALL ListCdagC_print(op%particles(flavor2),233)
  CALL ListCdagC_assign(op%list_swap, op%particles(flavor1)) !list_swap = particle
  op%particles(flavor1) = op%particles(flavor2)
  op%particles(flavor2) = op%list_swap
  !CALL ListCdagC_swap(op%particles(flavor1),op%particles(flavor2))
  !CALL ListCdagC_print(op%particles(flavor1),233)
  !CALL ListCdagC_print(op%particles(flavor2),233)

END SUBROUTINE ImpurityOperator_swap
!!***

!!****f* ABINIT/m_ImpurityOperator/ImpurityOperator_overlapIJ
!! NAME
!!  ImpurityOperator_overlapIJ
!!
!! FUNCTION
!!  Compute overlap between two 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=ImpurityOperator
!!  i=first flavor
!!  j=second flavor
!!
!! OUTPUT
!!  ImpurityOperator_overlapIJ=result
!!
!! SIDE EFFECTS
!!
!! NOTES
!!
!! PARENTS
!!  Will be filled automatically by the parent script
!!
!! CHILDREN
!!  Will be filled automatically by the parent script
!!
!! SOURCE

DOUBLE PRECISION FUNCTION ImpurityOperator_overlapIJ(op,i,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 'ImpurityOperator_overlapIJ'
!End of the abilint section

  TYPE(ImpurityOperator), INTENT(INOUT) :: op
  INTEGER               , INTENT(IN) :: i
  INTEGER               , INTENT(IN) :: j
!Local variables ------------------------------
!  TYPE(ListCdagC)       , POINTER    :: particle1 => NULL()
!  DOUBLE PRECISION, DIMENSION(:,:), POINTER :: list1 => NULL()
  INTEGER                            :: tail1
  DOUBLE PRECISION, DIMENSION(1:2)   :: CdagC_1
  INTEGER                            :: isegment

!  particle1 => op%particles(i) 
!  list1     => particle1%list
  tail1 = op%particles(i)%tail

  ImpurityOperator_overlapIJ = 0.d0
  IF ( tail1 .EQ. 0 .AND. op%particles(i)%list(0,C_) .EQ. 0.d0 ) THEN ! FULL
!    CALL CdagC_init(CdagC_1,0.d0,op%beta)
    CdagC_1(Cdag_) = 0.d0
    CdagC_1(C_   ) = op%beta

    ImpurityOperator_overlapIJ = ImpurityOperator_overlapSegFlav(op,CdagC_1,j)
  ELSE IF ( tail1 .GT. 0) THEN
    op%activeFlavor = i
    DO isegment = 1, tail1
        CdagC_1(:) = op%particles(i)%list(isegment,1:2)
        ImpurityOperator_overlapIJ = ImpurityOperator_overlapIJ &
                   + ImpurityOperator_overlapSegFlav(op,CdagC_1,j)
    END DO
  END IF

END FUNCTION ImpurityOperator_overlapIJ
!!***

!!****f* ABINIT/m_ImpurityOperator/ImpurityOperator_measDE
!! NAME
!!  ImpurityOperator_measDE
!!
!! FUNCTION
!!  measure double occupancy and interaction energy
!!
!! 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=ImpurityOperator
!!
!! OUTPUT
!!  DE=array accumulating duoble occupancy and energy
!!
!! SIDE EFFECTS
!!
!! NOTES
!!
!! PARENTS
!!  Will be filled automatically by the parent script
!!
!! CHILDREN
!!  Will be filled automatically by the parent script
!!
!! SOURCE

SUBROUTINE ImpurityOperator_measDE(op,DE)

!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 'ImpurityOperator_measDE'
!End of the abilint section

  TYPE(ImpurityOperator), INTENT(IN) :: op
  DOUBLE PRECISION, DIMENSION(:,:), INTENT(INOUT) :: DE
!Local variables ------------------------------
  DOUBLE PRECISION                              :: localD
  DOUBLE PRECISION                              :: totalE
  INTEGER                                       :: iflavor1
  INTEGER                                       :: iflavor2
  INTEGER                                       :: flavors

  IF ( .NOT. ALLOCATED(op%particles) ) &
    CALL ERROR("ImpurityOperator_measD : no particle set   ")

  totalE = 0.d0
  flavors = op%flavors
  DO iflavor1 = 1, flavors
    DO iflavor2 = iflavor1+1, flavors
      !localD = ImpurityOperator_overlapIJ(op,iflavor1,iflavor2) 
      localD = op%overlaps(iflavor2,iflavor1)
      DE(iflavor2,iflavor1) = DE(iflavor2,iflavor1) + localD  
      totalE = totalE + localD * op%mat_U(iflavor1,iflavor2)
    END DO
  END DO

  DE(1,1) = DE(1,1) + totalE

END SUBROUTINE ImpurityOperator_measDE
!!***

!!****f* ABINIT/m_ImpurityOperator/ImpurityOperator_cleanOverlaps
!! NAME
!!  ImpurityOperator_cleanOverlaps
!!
!! FUNCTION
!!  Compute from scratch all overlaps
!!
!! 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=ImpurityOperator
!!
!! OUTPUT
!!
!! SIDE EFFECTS
!!
!! NOTES
!!
!! PARENTS
!!  Will be filled automatically by the parent script
!!
!! CHILDREN
!!  Will be filled automatically by the parent script
!!
!! SOURCE

SUBROUTINE ImpurityOperator_cleanOverlaps(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 'ImpurityOperator_cleanOverlaps'
!End of the abilint section

  TYPE(ImpurityOperator), INTENT(INOUT) :: op
!Local variables ------------------------------
  INTEGER                                       :: iflavor1
  INTEGER                                       :: iflavor2
  INTEGER                                       :: flavors

  IF ( .NOT. ALLOCATED(op%particles) ) &
    CALL ERROR("ImpurityOperator_cleanOverlap : no particle set   ")

  flavors = op%flavors
  DO iflavor1 = 1, flavors
    DO iflavor2 = iflavor1+1, flavors
      op%overlaps(iflavor2,iflavor1) = ImpurityOperator_overlapIJ(op,iflavor1,iflavor2) 
    END DO
  END DO

END SUBROUTINE ImpurityOperator_cleanOverlaps
!!***

!!****f* ABINIT/m_ImpurityOperator/ImpurityOperator_measN
!! NAME
!!  ImpurityOperator_measN
!!
!! FUNCTION
!!  measure the number of electrons on flavor 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=ImpurityOperator
!!  flavor=the flavor
!!
!! OUTPUT
!!  ImpurityOperator_measN=number of electrons
!!
!! SIDE EFFECTS
!!
!! NOTES
!!
!! PARENTS
!!  Will be filled automatically by the parent script
!!
!! CHILDREN
!!  Will be filled automatically by the parent script
!!
!! SOURCE

DOUBLE PRECISION FUNCTION ImpurityOperator_measN(op,flavor)

!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 'ImpurityOperator_measN'
!End of the abilint section

  TYPE(ImpurityOperator), INTENT(IN) :: op
  INTEGER,      OPTIONAL, INTENT(IN) :: flavor
!Local variables ------------------------------
  DOUBLE PRECISION                   :: totalCdag
  DOUBLE PRECISION                   :: totalC
  INTEGER                            :: scanning
  INTEGER                            :: aF

  IF ( PRESENT(flavor) ) THEN
    aF = flavor
  ELSE
    aF = op%activeFlavor
  END IF

  IF ( aF .LE. 0 ) & 
    CALL ERROR("ImpurityOperator_measN : no active flavor     ")

  totalC    = (op%particles(aF)%list(0,Cdag_) - op%particles(aF)%list(0,C_) + op%beta) * .5d0 
  totalCdag = 0.d0

  DO scanning = 1, op%particles(aF)%tail
    totalCdag = totalCdag + op%particles(aF)%list(scanning,Cdag_)
    totalC    = totalC    + op%particles(aF)%list(scanning,C_   )
  END DO

  ImpurityOperator_measN = totalC - totalCdag

END FUNCTION ImpurityOperator_measN
!!***

!!****f* ABINIT/m_ImpurityOperator/ImpurityOperator_destroy
!! NAME
!!  ImpurityOperator_destroy
!!
!! FUNCTION
!!  destroy and deallocate
!!
!! 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=ImpurityOperator
!!
!! OUTPUT
!!
!! SIDE EFFECTS
!!
!! NOTES
!!
!! PARENTS
!!  Will be filled automatically by the parent script
!!
!! CHILDREN
!!  Will be filled automatically by the parent script
!!
!! SOURCE

SUBROUTINE ImpurityOperator_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 'ImpurityOperator_destroy'
!End of the abilint section

  TYPE(ImpurityOperator), INTENT(INOUT) :: op
!Local variables ------------------------------
  INTEGER                               :: IT

  IF ( ALLOCATED(op%particles) ) THEN
    DO IT = 1, op%flavors
      CALL ListCdagC_destroy(op%particles(IT))
    END DO
    DT_FREE(op%particles)
  ENDIF
  CALL ListCdagC_destroy(op%list_swap)
  FREEIF(op%mat_U)
  FREEIF(op%overlaps)
  FREEIF(op%updates)
  op%activeFlavor = 0
  op%flavors      = 0
  op%beta         = 0.d0
END SUBROUTINE ImpurityOperator_destroy
!!***

!!****f* ABINIT/m_ImpurityOperator/ImpurityOperator_getErrorOverlap
!! NAME
!!  ImpurityOperator_getErrorOverlap
!!
!! FUNCTION
!!  compute error on the overlap (numerical accumulation)
!!
!! 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=ImpurityOperator
!!
!! OUTPUT
!!  DE=save the error
!!
!! SIDE EFFECTS
!!
!! NOTES
!!
!! PARENTS
!!  Will be filled automatically by the parent script
!!
!! CHILDREN
!!  Will be filled automatically by the parent script
!!
!! SOURCE

SUBROUTINE ImpurityOperator_getErrorOverlap(op,DE)

!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 'ImpurityOperator_getErrorOverlap'
!End of the abilint section

  TYPE(ImpurityOperator), INTENT(INOUT) :: op
  DOUBLE PRECISION, DIMENSION(:,:), INTENT(INOUT) :: DE
!Local variables ------------------------------
  DOUBLE PRECISION                              :: localD1
  DOUBLE PRECISION                              :: localD2
  DOUBLE PRECISION                              :: totalE1
  DOUBLE PRECISION                              :: totalE2
  INTEGER                                       :: iflavor1
  INTEGER                                       :: iflavor2
  INTEGER                                       :: flavors

  IF ( .NOT. ALLOCATED(op%particles) ) &
    CALL ERROR("ImpurityOperator_getErrorOverlap : no particle set ")

  totalE1 = 0.d0
  totalE2 = 0.d0
  flavors = op%flavors
  DO iflavor1 = 1, flavors
    DO iflavor2 = iflavor1+1, flavors
      localD1 = ImpurityOperator_overlapIJ(op,iflavor1,iflavor2) 
      localD2 = op%overlaps(iflavor2,iflavor1)
      totalE1 = totalE1 + localD1 * op%mat_U(iflavor1,iflavor2)
      totalE2 = totalE2 + localD2 * op%mat_U(iflavor1,iflavor2)
    END DO
  END DO

  DE(2,2) = ABS(totalE1 - totalE2)

END SUBROUTINE ImpurityOperator_getErrorOverlap
!!***
!#ifdef CTQMC_CHECK
!!****f* ABINIT/m_ImpurityOperator/ImpurityOperator_doCheck
!! NAME
!!  ImpurityOperator_doCheck
!!
!! FUNCTION
!!  set the check mechanism
!!
!! 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=ImpurityOperator
!!  opt_check=1||3 do check
!!
!! OUTPUT
!!
!! SIDE EFFECTS
!!
!! NOTES
!!
!! PARENTS
!!  Will be filled automatically by the parent script
!!
!! CHILDREN
!!  Will be filled automatically by the parent script
!!
!! SOURCE

SUBROUTINE ImpurityOperator_doCheck(op,opt_check)

!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 'ImpurityOperator_doCheck'
!End of the abilint section

  TYPE(ImpurityOperator) , INTENT(INOUT) :: op
  INTEGER                , INTENT(IN   ) :: opt_check

  IF ( opt_check .EQ. 1 .OR. opt_check .EQ. 3 ) &
    op%doCheck = .TRUE.
END SUBROUTINE ImpurityOperator_doCheck
!!***

!!****f* ABINIT/m_ImpurityOperator/ImpurityOperator_checkOverlap
!! NAME
!!  ImpurityOperator_checkOverlap
!!
!! FUNCTION
!!  check the calculation of the overlap (very very slow routine)
!!  between Tmin and Tmax (c+ and c)
!!
!! 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=ImpurityOperator
!!  Tmin=c+
!!  Tmax=c
!!  iOverlap=input overlap (fast calculation)
!!  iflavor=active flavor
!!
!! OUTPUT
!!
!! SIDE EFFECTS
!!
!! NOTES
!!
!! PARENTS
!!  Will be filled automatically by the parent script
!!
!! CHILDREN
!!  Will be filled automatically by the parent script
!!
!! SOURCE

SUBROUTINE ImpurityOperator_checkOverlap(op, Tmin, Tmax, iOverlap, 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 'ImpurityOperator_checkOverlap'
!End of the abilint section

  TYPE(ImpurityOperator), INTENT(INOUT)  :: op
  DOUBLE PRECISION      , INTENT(IN   )  :: Tmin
  DOUBLE PRECISION      , INTENT(IN   )  :: Tmax
  DOUBLE PRECISION      , INTENT(IN   )  :: iOverlap
  INTEGER               , INTENT(IN   )  :: iflavor 
!Local variables ------------------------------
  INTEGER, PARAMETER                     :: size=10000000
  INTEGER                                :: imin
  INTEGER                                :: imax
  INTEGER                                :: imaxbeta
  INTEGER                                :: isegment
  INTEGER                                :: tail
  INTEGER(1), DIMENSION(1:size,1:2)      :: checktab 
  CHARACTER(LEN=4)                       :: a
  DOUBLE PRECISION                       :: dt
  DOUBLE PRECISION                       :: inv_dt
  DOUBLE PRECISION                       :: overlap
  DOUBLE PRECISION                       :: erreur
  DOUBLE PRECISION                       :: weight
  INTEGER :: try

  checktab = INT(0,1)
  overlap = 0.d0

  dt = op%beta / DBLE((size-1))
  inv_dt = 1.d0 / dt
  imin = INT(Tmin / dt + 0.5d0) + 1
  imax = INT(Tmax / dt + 0.5d0) + 1
  MODCYCLE(imax, size, imaxbeta)

  tail = op%particles(iflavor)%tail

  DO try = imin, MIN(imax,size)
    checktab(try,1)=INT(1,1)!IBSET(checktab(try,1),0)
  END DO

  IF ( imax .NE. imaxbeta ) THEN 
    DO try = 1, imaxbeta
      checktab(try,1)=INT(1,1)!IBSET(checktab(try,1),0)
    END DO
  END IF

  IF ( tail .NE. 0 ) THEN
    DO isegment=1, tail
      imin = INT(op%particles(iflavor)%list(isegment,Cdag_)* inv_dt + 0.5d0) + 1
      imax = INT(op%particles(iflavor)%list(isegment,C_   )* inv_dt + 0.5d0) + 1
      MODCYCLE(imax, size, imaxbeta)
      DO try = imin, MIN(imax,size)
        checktab(try,2)=INT(1,1)!IBSET(checktab(try,2),0)
      END DO
      IF ( imax .NE. imaxbeta ) THEN
        DO try = 1, imaxbeta
          checktab(try,2)=INT(1,1)!IBSET(checktab(try,2),0)
        END DO
      END IF
    END DO
  ELSE IF ( op%particles(iflavor)%list(0,C_) .EQ. 0.d0 ) THEN
    DO try = 1, size
      checktab(try,2)=INT(1,1)!IBSET(checktab(try,2),0)
    END DO
  END IF

  DO isegment = 1, size
    IF ( IAND(checktab(isegment,1),checktab(isegment,2)) .EQ. INT(1,1) ) &
      overlap = overlap + 1.d0
  END DO

  overlap = overlap * dt

  IF ( iOverlap .EQ. 0.d0 ) THEN
    erreur = ABS(overlap)
  ELSE
    erreur = ABS(overlap                - iOverlap)
  END IF
  weight = ABS(2.d0 * DBLE(tail) * dt - iOverlap)
  IF ( erreur .GT. weight  ) THEN 
    WRITE(a,'(I4)') INT(erreur*100.d0)
    CALL WARN("ImpurityOperator_checkOverlap : "//a//"%              ") 
  END IF
  IF ( iOverlap .LE. (2.d0 * DBLE(tail) * dt) ) &
    op%meanError = op%meanError + 1.d0
  op%checkNumber = op%checkNumber + 1.d0 !weight 

END SUBROUTINE ImpurityOperator_checkOverlap
!!***

!!****f* ABINIT/m_ImpurityOperator/ImpurityOperator_getError
!! NAME
!!  ImpurityOperator_getErro
!!
!! FUNCTION
!!  get error on computing the overlap
!!
!! 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=ImpurityOperator
!!
!! OUTPUT
!!  ImpurityOperator_getError=percentage error
!!
!! SIDE EFFECTS
!!
!! NOTES
!!
!! PARENTS
!!  Will be filled automatically by the parent script
!!
!! CHILDREN
!!  Will be filled automatically by the parent script
!!
!! SOURCE

DOUBLE PRECISION FUNCTION ImpurityOperator_getError(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 'ImpurityOperator_getError'
!End of the abilint section

  TYPE(ImpurityOperator), INTENT(IN) :: op
!Local variables ------------------------------
!  DOUBLE PRECISION :: tolerance
  DOUBLE PRECISION :: error

  IF ( op%doCheck .EQV. .TRUE. ) THEN
    error     = ABS(op%meanError/op%checkNumber) 
!  tolerance = ABS(op%tolerance/op%checkNumber) 
    ImpurityOperator_getError = error 
  ELSE
    ImpurityOperator_getError = 0.d0
  END IF
END FUNCTION ImpurityOperator_getError
!!***
!#endif

!!****f* ABINIT/m_ImpurityOperator/ImpurityOperator_printLatex
!! NAME
!!  ImpurityOperator_printLatex
!!
!! FUNCTION
!!  print in a latex format all the configuration
!!
!! 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=ImpurityOperator
!!  ostream=file stream
!!  isweep=current sweep number
!!
!! OUTPUT
!!
!! SIDE EFFECTS
!!
!! NOTES
!!
!! PARENTS
!!  Will be filled automatically by the parent script
!!
!! CHILDREN
!!  Will be filled automatically by the parent script
!!
!! SOURCE

SUBROUTINE ImpurityOperator_printLatex(op, ostream, isweep)

!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 'ImpurityOperator_printLatex'
!End of the abilint section

  TYPE(ImpurityOperator), INTENT(IN) :: op
  INTEGER               , INTENT(IN) :: ostream
  INTEGER               , INTENT(IN) :: isweep
!Local variables ------------------------------
  INTEGER                            :: flavors
  INTEGER                            :: iflavor
  INTEGER                            :: tail
  INTEGER                            :: it
  DOUBLE PRECISION                   :: C
  DOUBLE PRECISION                   :: Cdag
  INTEGER                            :: ordo
  INTEGER                            :: y
  INTEGER                            :: lines
  INTEGER                            :: letters
  DOUBLE PRECISION                   :: length

  flavors = op%flavors

  WRITE(ostream,'(A13)')    "\begin{frame}"
  WRITE(ostream,'(2x,A14)') "\begin{figure}"
  WRITE(ostream,'(4x,A28)') "\setlength{\unitlength}{1mm}"
  WRITE(ostream,'(4x,A23)') "\begin{picture}(104,90)"
  WRITE(ostream,'(6x,A29,I6,A2)') "\put(52,00){\makebox(0,0)[c]{",isweep,"}}"
  y = INT(90.d0/DBLE(flavors+1))
  DO iflavor = 1, flavors
    tail =  op%particles(iflavor)%tail
    ordo = iflavor * y
    lines = ordo - 1
    letters = ordo - 5
    WRITE(ostream,'(6x,A6,I2)') "%ligne", iflavor
    WRITE(ostream,'(6x,A41,I2,A16)') "\linethickness{0.5pt}\color{black}\put(2,",lines,"){\line(0,1){2}}"
    WRITE(ostream,'(6x,A7,I2,A24)')  "\put(2,",letters,"){\makebox(0,0)[c]{$0$}}"
    WRITE(ostream,'(6x,A7,I2,A18)')  "\put(2,",ordo,"){\line(1,0){100}}"
    WRITE(ostream,'(6x,A9,I2,A16)')  "\put(102,",lines,"){\line(0,1){2}}"
    WRITE(ostream,'(6x,A9,I2,A28)')  "\put(102,",letters,"){\makebox(0,0)[c]{$\beta$}}"
    DO it = 1, tail 
      Cdag = 2.d0+(op%particles(iflavor)%list(it,Cdag_)/op%beta*100.d0)
      C    = 2.d0+(op%particles(iflavor)%list(it,C_   )/op%beta*100.d0)
      length = C - Cdag
      IF ( op%particles(iflavor)%list(it,C_) .LE. op%beta ) THEN
        WRITE(ostream,'(8x,A9,I2)')             "%segments", it
        WRITE(ostream,'(8x,A37,F5.1,A1,I2,A13,F5.1,A2)') &
        "\linethickness{2pt}\color{black}\put(",Cdag,",",ordo,"){\line(1,0){",length,"}}"
        WRITE(ostream,'(8x,A5)')                "%Cdag"
        WRITE(ostream,'(8x,A12)')               "\color{blue}"
        WRITE(ostream,'(8x,A5,F5.1,A1,I2,A14)') "\put(",Cdag,",",ordo,"){\circle*{1}}"
        WRITE(ostream,'(8x,A2)')                "%C"
        WRITE(ostream,'(8x,A11)')               "\color{red}"
        WRITE(ostream,'(8x,A5,F5.1,A1,I2,A14)') "\put(",C,",",ordo,"){\circle*{1}}"
      ELSE
        WRITE(ostream,'(8x,A9,I2)')             "%segments", it
        WRITE(ostream,'(8x,A37,F5.1,A1,I2,A13,F5.1,A2)') &
        "\linethickness{2pt}\color{black}\put(",Cdag,",",ordo,"){\line(1,0){",102.d0-Cdag,"}}"
        WRITE(ostream,'(8x,A7,I2,A13,F5.1,A2)') "\put(2,",ordo,"){\line(1,0){",C-102.d0,"}}"
        WRITE(ostream,'(8x,A5)')                "%Cdag"
        WRITE(ostream,'(8x,A12)')               "\color{blue}"
        WRITE(ostream,'(8x,A5,F5.1,A1,I2,A14)') "\put(",Cdag,",",ordo,"){\circle*{1}}"
        WRITE(ostream,'(8x,A2)')                "%C"
        WRITE(ostream,'(8x,A11)')               "\color{red}"
        WRITE(ostream,'(8x,A5,F5.1,A1,I2,A14)') "\put(",C-100.d0,",",ordo,"){\circle*{1}}"
      END IF
    END DO
    IF ( tail .EQ. 0 .AND. op%particles(iflavor)%list(0,C_) .EQ. 0.d0 ) THEN 
      WRITE(ostream,'(8x,A9,I2)')      "%segments", it
      WRITE(ostream,'(8x,A39,I2,A18)') "\linethickness{2pt}\color{black}\put(2,",ordo,"){\line(1,0){100}}"
    END IF
  END DO
  WRITE(ostream,'(4x,A13)') "\end{picture}"
  WRITE(ostream,'(2x,A12)') "\end{figure}"
  WRITE(ostream,'(2x,A17)') "\transduration{0}"
  WRITE(ostream,'(A11)')    "\end{frame}"
END SUBROUTINE ImpurityOperator_printLatex
!!***

END MODULE m_ImpurityOperator
!!***
