A new three-layer method is presented that allows studying the catalytic effect of zeolites, and their selectivity, on chemical reaction. Three levels of approximation are combined. A small cluster described quantum mechanically represents the reactive center. The direct surrounding of this cluster is treated by molecular mechanics to explicitly consider the effects of chemical substitution on the reaction mechanism. Finally, the remaining of the crystal is modeled by the Madelung potential to take care of the long-range electrostatic interaction felt by the reactive part. The QM/MM method used is based on the previously developed LSCF/MM methodology originally created to study macromolecular systems. The electrostatic interaction between the Madelung potential and the QM/MM cluster is approximated by means of the expectation value of distributed multipole operators fitted to the electronic electrostatic potential. This method, after a very slight parameterization of the hybrid force field, is applied on the high-silica chabazite crystal (H-SSZ-13). The optimized geometrical structure is discussed with respect to experimental data and to other theoretical results. Finally, the relative proton affinity of the four oxygen atoms directly bonded to the aluminum atom substituting a silicon atom, is computed and the agreement with experimental evidences is quite satisfying.