Harmonic force field (FF) parameters for the active site of native azurin (AZ) have been developed using density functional theory (DFT)-based Becke's three-parameter hybrid exchange functional and the Lee-Yang-Parr correlation functional (B3LYP) method. The same computational protocol has also been applied to derive the FF parameters for the metal ion-substituted [Co(II) and Ni(II)] AZs. To validate the new set of FF parameters for the metal sites, molecular dynamics (MD) simulations on native, loop-contracted, and metal ion-substituted AZs have been carried out for 10 ns using AMBER parameters for the remaining part of the proteins. The average structure obtained from the MD simulation for native protein is akin to that of X-ray diffraction studies. Results from the in silico loop variation reveal that the active site of AZ is almost unaffected by the loop contraction in accordance with the previous experimental findings. However, the inherent hydrogen-bonded network of the metal site of AZ is affected by the loop contraction. Comparison of the average structures obtained from the MD simulations for the metal ion-substituted proteins with the corresponding X-ray diffraction structures shows that there are no major differences between two systems. Nevertheless, the metal ion binding site undergoes significant changes due to metal ion-substitution. Results clearly demonstrate the usefulness of a new set of FF parameters in the engineering and redesign of blue copper proteins.