Results are presented from a search for the Higgs boson decay H → Zγ, where Z → ℓ++ℓ−− with ℓ = e or μ. The search is performed using a sample of proton-proton (pp) collision data at a center-of-mass energy of 13 TeV, recorded by the CMS experiment at the LHC, corresponding to an integrated luminosity of 138 fb−1−1. Events are assigned to mutually exclusive categories, which exploit differences in both event topology and kinematics of distinct Higgs production mechanisms to enhance signal sensitivity. The signal strength μ, defined as the product of the cross section and the branching fraction \left[\sigma \left(\textrm{pp}\to \textrm{H}\right)\mathcal{B}\left(\textrm{H}\to \textrm{Z}\upgamma \right)\right] relative to the standard model prediction, is extracted from a simultaneous fit to the ℓ++ℓ−−γ invariant mass distributions in all categories and is measured to be μ = 2.4 ± 0.9 for a Higgs boson mass of 125.38 GeV. The statistical significance of the observed excess of events is 2.7 standard deviations. This measurement corresponds to \left[\sigma \left(\textrm{pp}\to \textrm{H}\right)\mathcal{B}\left(\textrm{H}\to \textrm{Z}\upgamma \right)\right]=0.21\pm 0.08 pb. The observed (expected) upper limit at 95% confidence level on μ is 4.1 (1.8), where the expected limit is calculated under the background-only hypothesis. The ratio of branching fractions \mathcal{B}\left(\textrm{H}\to \textrm{Z}\upgamma \right)/\mathcal{B}\left(\textrm{H}\to \upgamma \upgamma \right) is measured to be 1.5−0.6+0.71.5−0.6+0.7, which agrees with the standard model prediction of 0.69 ± 0.04 at the 1.5 standard deviation level.