anchor.diff:  !    call print_Mat(Dens(:,:,iSpin),dim_1e,5,"DensX")
     enddo
   
  +  ! Damping    
  +  if(DoDamp)then
  +    Dens = Damp*Dens + (1.0d0-Damp)*Densold
  +!    call calc_DynDamp()
  +  elseif(DynDamp)then
  +    call calc_DynDamp()
  +    Dens = Damp*Dens + (1.0d0-Damp)*Densold
  +  endif
  +
     Drmsd = sqrt(Drmsd/dble(dim_1e*dim_1e*Spins))
   
  -
   end subroutine calc_Dens
   
  +
  +!#############################################
  +!#       Calculate Dynamic Damping Factor
  +!#############################################
  +subroutine calc_DynDamp()
  +  use module_data, only      : Fock,Coef,S12,dim_1e,Spins,Occ
  +  use module_data, only      : Dens,Densold,Hcore,Fockold,Damp
  +  implicit none
  +  double precision          :: E1a, E2a, E1b, E2b, E2ab, temp,tempa,tempb
  +  integer                   :: i,j,iSpin
  +
  +  E1a  = 0.0d0
  +  E1b  = 0.0d0
  +  E2a  = 0.0d0
  +  E2b  = 0.0d0
  +  E2ab = 0.0d0
  +
  +!        Eelec = Eelec + Dens(i,j,iSpin)                &
  +!                      *(Hcore(i,j) + Fock(i,j,iSpin))  &
  +!                      *(1.0d0/dble(Spins))
  +
  +  do iSpin=1,Spins
  +    do i=1,dim_1e
  +      do j=1,dim_1e
  +        E1a   = E1a + Dens(i,j,iSpin)*(Hcore(i,j)*(1.0d0/dble(Spins)))
  +        E2a   = E2a + Dens(i,j,iSpin)*((Fock(i,j,iSpin))*(1.0d0/dble(Spins)))
  +
  +        E1b   = E1b + Densold(i,j,iSpin)*(Hcore(i,j)*(1.0d0/dble(Spins)))
  +        E2b   = E2b + Densold(i,j,iSpin)*((Fockold(i,j,iSpin))*(1.0d0/dble(Spins)))
  +
  +        E2ab  = E2ab + Dens(i,j,iSpin)*((Fockold(i,j,iSpin))*(1.0d0/dble(Spins)))
  +
  +      enddo
  +    enddo
  +  enddo
  +
  +
  +  tempa = E1a + E2a
  +  tempb = E1b + E2b
  +  temp  = (E1a-E1b+E2ab-2.0d0*E2b)/(4.0d0*E2ab-2.0d0*E2a-2.0d0*E2b)
  +  damp = min(abs(temp),1.0d0)
  +
  +!  write(*,*) E1a,"-",E1b,"+",E2ab,"-2",E2b
  +!  write(*,*) "4",E2ab,"-2",E2a,"-2",E2b
  +!  write(*,*) tempa,tempb,E2ab,temp 
  +
  +end subroutine calc_DynDamp
   
   
   !#############################################
  
positive.code:     Fock(1:dim_1e,1:dim_1e,iSpin) = Hcore(1:dim_1e,1:dim_1e)
    enddo
  
    do i=0,dim_1e-1
      do j=0,dim_1e-1
        do k=0,dim_1e-1
  
negative.code: 
      enddo
    enddo
  !$OMP END DO
  !$OMP END PARALLEL
  
    ! Damping    
  !  if(DoDamp)then
  !    Fock = Damp*Fock + (1.0d0-Damp)*Fockold
  !  endif
  
  !  do iSpin=1,Spins
  !    call print_Mat(Fock(:,:,iSpin),dim_1e,4,"Fock")
  !  enddo
  
  !  deallocate(Fockold)
  
  end subroutine calc_Fock
  
  