Molecular optical susceptibilities are calculated by deriving equations of,motion for the single electron reduced density matrix, and solving them using the time dependent Hartree-Fock (TDHF) approximation. The present approach focuses directly on the dynamics of the charges in real space and completely avoids the tedious summations over molecular eigenstates. It further maps the system onto a set of coupled harmonic oscillators. The density matrix clearly shows the electronic structures induced by the external field, and how they contribute to the optical response. The-method is applied to calculating;the frequency-dispersed optical susceptibility chi((3)) of conjugated linear polyenes, starting with the Pariser-Parr-Pople (PPP) model. Charge density wave (CDW) like fluctuations and soliton pair like local bond-order fluctuations are shown to play important roles in the optical response of these systems.