Computational protein design (CPD) is an important tool for biotechnology still under development. Early applications led to proteins with novel ligand-binding functions, novel enzyme activity, and proteins that were completely "redesigned": around 2/3 of their sequence was mutated, yet their structure and stability were retained. In the last few years, CPD has allowed the creation of new protein folds, completely new enzymes, and the assembly or deassembly of multiprotein complexes. CPD methods are mainly characterized by (a) the energy function, (b) the description of the folded protein's conformational space, (c) the treatment of the unfolded state, and (d) the search method used to explore sequences and conformations. Graphical model and cost function network have been recently introduced as new search methods for CPD search allowing to find the global minimum energy conformation (with optimality proof). However, it can solved problems with less than 1003 and 484 mutations respectively with Rosetta and Amber forcefield and Dunbrack and Tuffery rotamers library. In this work, we introduce a new search method for CPD dedicated to large instances. The method, based on Variable Neighborhood Search (VNS), uses (partial) tree search in order exhibit the solutions with the lowest energy. In order to improve its behavior, we implemented a new heuristic taking advantage of the euclidean information provided by the pdb structure for neighbor variable selections. This last one used in conjunction with logarithmic size incrementation of the neighborhood size improves significantly the VNS behaviors. The resulting algorithm, called euclidean VNS, is more robust. It also outperforms Replica-Exchange Method for global minimum energy search. Furthermore, the method can provide direct correlation between energy improvement and protein structure, allowing to identify energetic hot spots in the protein backbone.