We report ligand field molecular mechanics, density functional theory and semi-empirical studies on the binding of Cu(II) to GlyHisLys (GHK) peptide. Following exhaustive conformational searching using molecular mechanics, we show that relative energy and geometry of conformations are in good agreement between GFN2-xTB semi-empirical and B3LYP-D DFT levels. Conventional molecular dynamics simulation of Cu-GHK shows the stability of the copper-peptide binding over 100ps trajectory. Four equatorial bonds in 3N1O coordination remain stable throughout simulation, while a fifth in apical position from Cterminal carboxylate is more fluxional. We also show that the automated conformer and rotamer search algorithm CREST is able to correctly predict the metal binding position from a starting point consisting of separated peptide, copper and water.