The implementation of the analytical energy gradient method for the calculation of first order one electron properties using closed shell configuration interaction wave functions of the single reference plus single and double substitutions (CI-SD) type is discussed. The method used is based on the Z vector formalism of Handy and Schaefer which is readily extended to allow a full orbital optimization to be carried out for a given CI-SD wave function. The results of comparative test calculations are reported for the HF, H 2O, CO, HCN, and O3 molecules, for which the dipole and quadrupole moments and the electric field gradients at the nuclei have been calculated by both the expectation value and the energy derivative formalisms using several standard basis sets. The effects of orbital optimization on the above properties at the equilibrium geometries as well as at a range of distorted geometries for HF, are also discussed. It is found that agreement between the different formalisms is best when the reference state in the CI expansion is strongly dominant. The DHS scheme of Pulay has been incorporated into the orbital optimization method and has been found to be efficient in generating the fully optimized CI-SD/MCSCF wave functions.