Let \(\Omega \subseteq \mathbb {R}^{d}\) be open, A a complex uniformly strictly accretive \(d\times d\) matrix-valued function on \(\Omega \) with \(L^{\infty }\) coefficients, b and c two d-dimensional vector-valued functions on \(\Omega \) with \(L^{\infty }\) coefficients and V a locally integrable nonnegative function on \(\Omega \) . Consider the operator \({{\mathscr {L}}}^{A,b,c,V}=-\textrm{div}\,(A\nabla \cdot ) + \left\langle \nabla , \overline{b}\right\rangle - \textrm{div}\,(c \, \cdot ) + V \) with mixed boundary conditions on \(\Omega \) . We extend the bilinear inequality that Carbonaro and Dragičević proved in [Bilinear embedding for Schrödinger-type operators with complex coefficients. Publ. Mat. (to appear)] in the special cases when \(b=c = 0\) , previously proved in (Calc Var Part Differ Equ 59(3):36, Paper No. 104, 2020) when \(V=0\) as well. As a consequence, we obtain that the solution to the parabolic problem \(u^{\prime }(t)+{{\mathscr {L}}}^{A,b,c,V}u(t)=f(t)\) , \(u(0)=0\) , has maximal regularity in \(L^{p}(\Omega )\) , for all \(p>1\) such that A satisfies the p-ellipticity condition that Carbonaro and Dragičević introduced in (J Eur Math Soc 22(10):3175–3221, 2020) and b, c, V satisfy another condition that we introduce in this paper. Roughly speaking, V has to be “big” with respect to b and c. We do not impose any conditions on \(\Omega \) , in particular, we do not assume any regularity of \(\partial \Omega \) , nor the existence of a Sobolev embedding.