The rare B(s) --> gammagamma decay has been investigated in the framework of the two-Higgs-doublet model. It is shown that the branching ratio has a great enhancement in comparison to the standard model (SM) prediction. We have found that the ratio of the square of the CP = 1 to CP = -1 amplitudes, \A+\2/\A-\2, which is always less than one in the SM, becomes greater than one for some values of M(H)/M(t). We have also calculated the direct CP-violating parameter epsilon'(gammagamma(-)). It is found that for K --> gammagamma decay its value is greater than the one which follows from the SM, while for B(s) --> gammagamma decay both models have the same predictions. Finally, we compare our results on B(s) --> gammagamma decay with the results of B(s) --> mu+mu- and B(s) --> mu+mu- decays which are of the same order of the coupling constant.