We present ab initio calculations of frequency-dependent linear and nonlinear optical responses based on real-time time-dependent density functional theory for arbitrary photonic molecules. This approach is based on an extension of an approach previously implemented for a linear response using the electronic structure program SIESTA. Instead of calculating excited quantum states, which can be a bottleneck in frequency-space calculations, the response of large molecular systems to time-varying electric fields is calculated in real time. This method is based on the finite field approach generalized to the dynamic case. To speed the nonlinear calculations, our approach uses Gaussian enveloped quasimonochromatic external fields. We thereby obtain the frequency-dependent second harmonic generation beta(-2omega;omega,omega), the dc nonlinear rectification beta(0;-omega,omega), and the electro-optic effect beta(-omega;omega,0). The method is applied to nanoscale photonic nonlinear optical molecules, including p-nitroaniline and the FTC chromophore, i.e., 2-[3-Cyano-4-(2-{5-[2-(4-diethylamino-phenyl)-vinyl]-thiophen-2-yl}-vinyl)-5,5-dimethyl-5H-furan-2-ylidene]-malononitrile, and yields results in good agreement with experiment.