Optimization in the normal mode coordinates has been established as a useful tool for modeling of vibrational spectra ( J. Chem. Phys. 2002, 117, 4126). In this work the algorithm is extended with the aid of harmonic penalty functions to allow for multiple restraints of geometry parameters, such as bond lengths, bond and dihedral angles, and for simultaneous optimization of more molecules. Additionally, geometry optimization when atomic nuclei are maintained on the constant electrostatic potential surface was implemented and its applications for solvent models are discussed. Model systems include small molecules, water cluster, antiparallel beta-sheet peptide containing intermolecular hydrogen bonds, periodic alpha-helix and a parallel beta-sheet segments. The normal mode method provided better numerical stability than the conventional redundant internal coordinates, especially for weakly hydrogen-bonded systems, while the speed of the optimization was found similar as for the Cartesian coordinates.