We present a detailed theoretical study of nonfactorizable contributions of the charm-quark loop to the amplitude of the $B_s \to \gamma \gamma$ decay. This contribution involves the B-meson three-particle Bethe-Salpeter amplitude, for which we take into account constraints from analyticity and continuity. The charming-loop contribution of interest may be described as a correction to the Wilson coeffcient $C_{7\gamma}$, $C_{7\gamma} \to C_{7\gamma} (1 + \delta C_{7\gamma})$. We calculate an explicit dependence of $\delta C_{7\gamma}$ on the parameter $\lambda_{B_s}$. Taking into account all theoretical uncertainties, $\delta C_{7\gamma}$ may be predicted with better than $10 \%$ accuracy for any given value of $\lambda_{B_s}$. For our benchmark point $\lambda_{B_s} = 0.45$ GeV, we obtain $\delta C_{7\gamma} = 0.045 \pm 0.004$. Presently, $\lambda_{B_s}$ is not known with high accuracy, but its value is expected to lie in the range $0.3 \leq \lambda_{B_s} (GeV) \leq 0.6$. The corresponding range of $\delta C_{7\gamma}$ is found to be $0.02 \leq \delta C_{7\gamma} \leq 0.1$. One therefore expects the correction given by charming loops at the level of at least a few percent.