  By employing the QCD factorization approach, we calculated the
next-to-leading order new physics contributions to the branching ratios, CP
asymmetries, isospin and U-spin symmetry breaking of the exclusive decays $B
\to V \gamma$ ($V=K^*, \rho$), induced by the charged Higgs penguins in the
general two-Higgs-doublet models. Within the considered parameter space, we
found that (a) the new physics corrections to the observables are generally
small in the model I and model III-A, moderate in model II, but large in model
III-B; (b) from the well measured branching ratios and upper limits, a lower
bound of $\mhp > 200$ GeV in model II was obtained, while the allowed range of
$\mhp$ in model III-B is $ 226 \leq \mhp \leq 293 $ GeV; these bounds are
comparable with those from the inclusive $B \to X_s \gamma$ decay; (c)the NLO
Wilson coefficient $C_7(m_b)$ in model III-B is positive and disfavered by the
measured value of isospin symmetry breaking $\Delta_{0-}^{exp} (K^*\gamma) =
(3.9 \pm 4.8)%$, but still can not be excluded if we take the large errors into
account; (d) the CP asymmetry $\acp(B \to \rho \gamma)$ in model III-B has an
opposite sign with the one in the standard model (SM), which may be used as a
good observable to distinguish the SM from model III-B; (e) the isospin
symmetry breaking $\Delta(\rho\gamma)$ is less than 10% in the region of
$\gamma = [ 40 \sim 70]^\circ$ preferred by the global fit result, but can be
as large as 20 to 40% in the regions of $\gamma \leq 10^\circ$ and $\gamma \geq
120^\circ$. The SM and model III-B predictions for $\Delta(\rho \gamma)$ are
opposite in sign for small or large values of the CKM angles; (f) the U-spin
symmetry breaking $\Delta U(K^*,\rho)$ in the SM and the general
two-Higgs-doublet models is generally small in size: $\sim 10^{-7}$.
