A previous article proposed an electronic structure-based polarizable potential, called the explicit polarization (X-POL) potential, to treat many-body polarization and charge delocalization effects in polypeptides. Here, we present a variational version of the X-POL potential, in which the wave function of the entire molecular system is variationally optimized to yield the minimum total electronic energy. This allows the calculation of analytic gradients, a necessity for efficient molecular dynamics simulations. In this paper, the detailed derivations of the Fock matrix and analytic force are presented and discussed. The calculations involve a double self-consistent-field procedure in which the wave function of each fragment is self-consistently optimized in the presence of other fragments, and in addition the polarization of the entire system is self-consistently optimized. The variational X-POL potential has been implemented in the Chemistry at Harvard Molecular Mechanics (CHARMM) package and tested successfully for small model compounds.