anchor.diff:      max_e_pert = 0d0
       
       do istate=1,N_states
  -      val = vect(istate, p1)
  +      val = vect(istate, p1) + vect(istate, p1)
         delta_E = E0(istate) - Hii
  +      tmp = dsqrt(delta_E * delta_E + val * val)
         if (delta_E < 0.d0) then
  -        e_pert = 0.5d0 * (-dsqrt(delta_E * delta_E + 4.d0 * val * val) - delta_E)
  -      else
  -        e_pert = 0.5d0 * ( dsqrt(delta_E * delta_E + 4.d0 * val * val) - delta_E)
  +        tmp = -tmp
         endif
  +      e_pert = 0.5d0 * ( tmp - delta_E)
         pt2(istate) += e_pert
         if(dabs(e_pert) > dabs(max_e_pert)) max_e_pert = e_pert
       end do
  
positive.code:   
    integer, intent(in) :: i,j,k,l
    double precision, external :: get_mo_bielec_integral
    
    integral8 = get_mo_bielec_integral(i,j,k,l,mo_integrals_map)
  end function
  
  
  
negative.code:       
        do istate=1,N_states
          delta_E = E0(istate) - Hii
          val = mat(istate, p1, p2) + mat(istate, p1, p2) 
          tmp = dsqrt(delta_E * delta_E + val * val)
          if (delta_E < 0.d0) then
            tmp = -tmp
          endif
          e_pert = 0.5d0 * ( tmp - delta_E)
          pt2(istate) = pt2(istate) + e_pert
          max_e_pert = min(e_pert,max_e_pert)
        end do
        
        if(dabs(max_e_pert) > buf%mini) then
  