The recently developed Variable Neighborhood Search (VNS) meta-heuristic for combinatorial and global optimization is outlined together with its specialization to the problem of finding extremal graphs with respect to one or more invariants, and the corresponding program (AGX). We illustrate the potential of the VNS algorithm on the example of energy E, a graph invariant which (in the case of molecular graphs of conjugated hydrocarbons) corresponds to the total -electron energy. Novel lower and upper bounds for E are suggested by AGX and several conjectures concerning (molecular) graphs with extremal E-values put forward. Moreover, most of the bounds are proved to hold.