The method of calculating distributed polarizabilities is extended to the first and second dipole hyperpolarizabilities, in order to describe more accurately the molecular response to strong and inhomogeneous external time-dependent electric fields. The dipolar response is expressed in terms of both potential related charge-density response functions and electric field related dipole-density response functions. The macroscopic linear, quadratic, and cubic optical dipole susceptibilities of molecular crystals are expressed in terms of the distributed (hyper) polarizabilities. This formulation differs from previous theories using distributed dipoles in that it allows for a rigorous treatment of both local induced dipoles and charge flow between different regions of the molecule. As an example, the distributed polarizabilities and first hyperpolarizabilities of urea at the self-consistent-field level are used to calculate the linear and quadratic susceptibilities of the urea crystal. The linear susceptibility does not differ substantially from that calculated with previous less rigorous models for distributed response, but the quadratic susceptibility is about 50% of that calculated with previous models. This indicates that the present treatment of distributed response should give a quadratic susceptibility in good agreement with experimental data, once the effects of electronic correlation, frequency dispersion, and the permanent crystal field are taken into account. (C) 2000 American Institute of Physics. [S0021-9606(00)30111-8].