Multi-wavelength spectral energy distributions of Low Mass X-ray Binaries in the hard state are determined by the emission from a jet, for frequencies up to mid-infrared, and emission from the accretion flow in the optical to X-ray range. In the last years, the flat radio-to-mid-IR spectra of Black Hole (BH) X-ray binaries was described using the internal shocks model, which assumes that the fluctuations in the velocity of the ejecta along the jet are driven by the fluctuations in the accretion flow, described by the X-ray Power Density Spectrum (PDS). In this work we attempt to apply this model for the first time to a Neutron Star (NS) LMXB, i.e. 4U 0614+091. We used the multi-wavelength data set obtained in 2006, comprising data from radio to X-ray, and applied a model which includes an irradiated disc model for the accretion flow and an updated version of the internal shocks code for the ejection. The new version of the code allows to change the geometry of the jet for the case of non-conical jets. Only two alternative scenarios provide a satisfactory description of the data: using the X-ray PDS but in a non-conical geometry for the jet, or either using a conical geometry but with a flicker-noise PDS. Both scenarios would imply some differences with the results obtained with similar models on BH X-ray binaries, shedding light on the possibility that jets in NS and BH binaries might somehow have a different geometry or a different coupling with the accretion flow.