In my PhD dissertation I report the applications of the Rigorous Local Field Theory and its evolutionary extensions for the description of linear and second-order nonlinear optical properties of organic crystals. The methodology capabilities were thoroughly investigated basing mainly on the molecular crystal of MNA. The in-depth considerations concerning the chromophore geometry, effects of the in-crystal electric field, computational methodology and ZPVA correction were discussed in a systematic way. This allowed to modernize the RLFT approach applications, scattered so far throughout the literature, and the tests performed for a vast family of nonlinear crystals proved the high predictive power of the methodology. The first step was the inclusion of the polarizing electric field. The uniform dipole field approximation was proved to give satisfactory agreement with experimental macroscopic electric properties. The enhancement of the molecular properties occur following the calculated dressing electric field. An improvement was obtained by the use of point charges field introducing non-uniformity which is favorable for larger molecules. The use of the neutral spherical surrounding showed, especially in the case of MNA, that the charges lying in the outermost regions of the shell introduce an erroneous condenser-like electric field. The neglect of the demand of surrounding electroneutrality leads to problems with convergence which is expected from the viewpoint of divergent manner of the dipole-dipole, charge-dipole and charge-charge interactions in an infinite perfect crystal lattice. The usefulness of this approach was tested/proved by performing calculations for total of 10 selected organic crystals. Several solutions for the optimum chromophore geometry, which is of central importance for the RLFT methodology, were considered. The non-optimized X-ray geometries usually give too small electric response as the protons are not localized precisely enough and perturb the precise localization of the other nuclei. The quality of geometries obtained in diffraction experiments may vary to a large extent depending on many conditions, ie. the crystal quality, measurement setup, data refinement strategy etc. The isolated molecule geometry optimization neglects the in-crystal surrounding of chromophore and sometimes leads to drastic conformational changes as evidenced in case of DAN. The periodic boundary conditions optimization is also constraint to limitations of the functional used. In the case of B3LYP used in the course of the study a systematic overestimation of the bond length alternation was observed which leads, for this class of compounds, to slight overestimations of the electric response. Further methodology improvements towards improvement of the X-ray geometry is highly desirable as the cost of experiment is relatively low and simple when compared ie. to neutron diffraction. One of the options might be optimization in the point charges field and has to be explored more thoroughly in the future. The RLFT methodology allowed also to assess the performance/influence of the isolated chromophore computational method choice on the macroscopic electric properties. The electron correlation effect was assessed for the MNA molecular crystal. The MP2 method was proved to recover a large part of the electron correlation. The performance of connection of the MP2 method for the static electric properties with B3LYP for frequency dispersion was comparable to the most sophisticated method used - CCSD. The performance of different XC functionals allowed to conclude that their performance correlates strongly with the amount of the HF exchange. This is a consequence of the experimental excitation energies reproduction. The lowest responses were obtained for pure HF method which is due to overestimating of excitation energies. The completely opposite side of the scale was occupied by the BLYP functional (with no HF exchange) and giving the strongest overestimations of the molecular responses. The range corrected methods, LC-BLYP and CAM-B3LYP, yielded rather underestimations however their usage is advisable for molecular systems of higher spacial n-conjugated backbone extent. All the tested functionals were shown to be much less sensitive towards the geometry changes than the wavefunction methods MP2, MP4 and CCSD. In the case of small NLO chromophores studied during my project the usefulness of B3LYP method for the frequency dispersion was supported to be excellent. The magnitude of the ZPVA correction was assessed for MNA with B3LYP method. It was shown to attain as much as 5% and 10-15% for dynamic x (1) and x (2) tensor components, respectively. These values may be treated as a rule of thumb, however, an automated procedure of calculations will extend the range of compounds treated with the finite distortion approach and lead to more general conclusions. The ionic organic crystals of considerable NLO activity were the subject of calculations in this study for the first time in the literature. The reduction of the microscopic response by the in-crystal dressing field was found in contrast with all the molecular crystals investigated in this thesis. The experimentally determined sequence with respect to non-resonant xm values: DAPSH>DSTMS~DAST was reproduced with very good accuracy. To sum up, the investigations conducted throughout this PhD project allowed to review and improve the Rigorous Local Field Theory approach. The current state-of-the-art allows to apply this multi-scale methodology for the interpretative and predictive purposes with high degree of confidence. A combined effort of experimentalists and theoretical chemists will contribute to further developments of this approach. These should include the modification of the partitioning scheme by using e.g. Quantum Theory of Atoms in Molecules (QTAIM) due to Bader and coworkers, vibrational contributions to the optical rectification process and description of third-order nonlinearities as well as the optical activity of crystals.