fix mismatch between intra and unbound energy

This bug did not affect poses but it affected reported energies.
In Vina (not in AD4), the UNBOUND is INTRA of the first pose.
 - macrocycle glue potential was added to UNBOUND but not to
   INTRA in the output pdbqt (during the search it added to
   the cost function as it should be). This would cause
   discrepancies of about 10 kcal/mol per glue pair.
 - poses were not resorted after refining with non_cache, allowing
   the UNBOUND energy to be the INTRA from a pose other than the first.
   This affected many poses by a few tenths of a kcal/mol, and some
   poses by a few kcal/mol.
 - After refinement with non_cache, it was possible to have movable
   atoms just outside the grid box because of a positive tolerance
   margin set in non_cache.h. The out-of-box penalty used during
   refinement is 100/angstrom in the first refinement iteration,
   but the unbound (intramolecular) energy is calculated with
   1,000,000/angstrom. Thus, in rare cases, flex sidechain atoms
   would be refined to be just outside the grid box with a tiny
   out-of-box penalty, then sorted using the tiny penalty,
   and finally rescored with an UNBOUND energy using the large
   penalty, leading to UNBOUND that can reach 100 kcal/mol.

src/lib/model.cpp

@@ -875,15 +875,6 @@ fl model::eval_intramolecular(const precalculate_byatom& p, const igrid& ig, con
 		}
 	}
 
-	// glue_i - glue_i and glue_i - glue_j interactions (no cutoff)
-	VINA_FOR_IN(i, glue_pairs) {
-		const interacting_pair& pair = glue_pairs[i];
-		fl r2 = vec_distance_sqr(coords[pair.a], coords[pair.b]);
-		fl this_e = p.eval_fast(pair.a, pair.b, r2);
-		curl(this_e, v[2]);
-		e += this_e;
-	}
-
 	return e;
 }
 

src/lib/non_cache.h

@@ -33,7 +33,7 @@ struct non_cache : public igrid {
 	virtual fl eval      (const model& m, fl v) const; // needs m.coords // clean up
 	virtual fl eval_intra(      model& m, fl v) const;
 	virtual fl eval_deriv(      model& m, fl v) const; // needs m.coords, sets m.minus_forces // clean up
-	bool within(const model& m, fl margin = 0.0001) const;
+	bool within(const model& m, fl margin = -0.0001) const; // negative margin ensures we stay inside (thanks Andreas!)
 	fl slope;
 private:
 	szv_grid sgrid;

src/lib/vina.cpp

@@ -944,6 +944,7 @@ void Vina::global_search(const int exhaustiveness, const int n_poses, const doub
 				}
 			}
 
+			poses.sort(); // order often changes after non_cache refinement
 			m_model.set(poses[0].c);
 			if (m_no_refine || !m_receptor_initialized)
 				intramolecular_energy = m_model.eval_intramolecular(m_precalculated_byatom, m_grid, authentic_v);
@@ -973,8 +974,10 @@ void Vina::global_search(const int exhaustiveness, const int n_poses, const doub
 			}
 		}
 
-		// Since pose.e contains the final energy, we have to sort them again
-		poses.sort();
+		// In AD4, the unbound energy is intra for each pose, so order may have changed
+		// The order does not change in Vina because unbound is intra of the 1st pose
+		if (m_sf_choice == SF_AD42)
+			poses.sort();
 
 		// Now compute RMSD from the best model
 		// Necessary to do it in two pass for AD4 scoring function