anchor.diff:    !$OMP END PARALLEL
     call end_parallel_job(zmq_to_qp_run_socket, 'pt2')
     double precision :: E0, avg, eqt
  -  double precision :: sumabove(comb_teeth), sum2above(comb_teeth), Nabove(comb_teeth)
     call do_carlo(tbc, Ncomb, comb, pt2_detail, sumabove, sum2above, Nabove)
  -  tbc(0) = 0
     !END LOOP?
     integer :: tooth
     !-8.091550677158776E-003
     call get_first_tooth(computed, tooth)
  -  E0 = sum(pt2_detail(1,:first_det_of_teeth(tooth)-1))
  -  avg = E0 + (sumabove(tooth) / Nabove(tooth))
  -  eqt = sqrt(1d0 / (Nabove(tooth)-1) * abs(sum2above(tooth) / Nabove(tooth) - (sumabove(tooth)/Nabove(tooth))**2))
  -  print *, "PT2 ", avg, "+/-", eqt
  +  print *, "TOOTH ", tooth 
  +  
  +  !!! ASSERT
  +  do i=1,first_det_of_teeth(tooth)-1
  +    if(not(computed(i))) stop "deter non calc"
  +  end do
  +  logical :: ok
  +  ok = .false.
  +  do i=first_det_of_teeth(tooth), first_det_of_teeth(tooth+1)-1
  +    if(not(computed(i))) ok = .true.
  +  end do
  +  if(not(ok)) stop "not OK..."
  +  !!!!!
  +
  +  if(Nabove(tooth) >= 30) then
  +    E0 = sum(pt2_detail(1,:first_det_of_teeth(tooth)-1))
  +    avg = E0 + (sumabove(tooth) / Nabove(tooth))
  +    eqt = sqrt(1d0 / (Nabove(tooth)-1) * abs(sum2above(tooth) / Nabove(tooth) - (sumabove(tooth)/Nabove(tooth))**2))
  +    time = omp_get_wtime()
  +    print "(A, 5(E15.7))", "PT2stoch ", time - time0, avg, eqt, Nabove(tooth)
  +  else
  +    print *, Nabove(tooth), "< 30 combs"
  +  end if
  +  tbc(0) = 0
  +  end do
  +
     pt2 = 0d0
   end subroutine
   
  
positive.code:     first_det_of_teeth(i) = pt2_find(stato, cweight)
      stato -= comb_step
    end do
  print *, first_det_of_teeth(1), first_det_of_comb
    if(first_det_of_teeth(1) /= first_det_of_comb) stop "comb provider"
  END_PROVIDER
  
  
  
negative.code:     myVal = 0d0
      myVal2 = 0d0
      do j=comb_teeth,1,-1
        !if(pt2_detail(1, dets(j)) == -1d0) print *, "uncalculatedidified", dets(j), pt2_detail(1, dets(j)-1:dets(j)+1)
        myVal += pt2_detail(1, dets(j)) / weight(dets(j)) * comb_step
        sumabove(j) += myVal
        sum2above(j) += myVal**2
  