Dear Patrizio,
I am glad that I understood the problem and that you were able to solve it with some of my suggestions. Note that my third option also included trying a different vdW correction. I was thinking especially about the family of nonlocal vdW-DF functionals. For those, the dispersion correction is baked into the functional, so it is calculated at every electronic step and not at the end of the calculation, so if it runs out of memory, you would at least notice this much earlier. However, I completely understand that you might need to stick to the DFT-D3 correction method for many reasons.
Now to the ENCUT question:
First of all, you are completely right that the quantity of interest (in your case, the phonon frequencies) needs to be converged with respect to ENCUT. However, 1000 eV seems very high for the chosen pseudopotentials. While it is correct that convergence studies might result in a necessary ENCUT > ENMAX, a > 100% increase might lead to problems since the pseudopotentials were not designed and tested for such a large increase. Results are probably fine, but you might run into many different issues, including convergence problems.
For phonon calculations, a lot of different parameters need to be converged, other than ENCUT. Especially convergence of the supercell size is important. You can see some examples in our phonon tutorials. If I read your POSCAR correctly, this is a 2x2x2 supercell? That might be a bigger problem regarding imaginary modes than the lower ENCUT value. But you are right in testing both. Since your unit cell is quite large, a further increase might be computationally challenging, but if you are still getting imaginary frequencies, it could be necessary. (Note that imaginary frequencies could also point to a dynamically unstable crystal structure and a possible phase transition. I have no experience with molecular systems, however, so I am hesitant to give advice here.)
Where you are a bit mistaken, if I read your last post correctly, is converging the total energy with respect to ENCUT for the primitive cell. Total energies are meaningless in DFT; only energy differences are important. For example, the energy difference between two volumes or between the unperturbed system and the one with a phonon displacement. Trying to converge the total energy for a single system is thus generally a bad idea, regardless of what parameter you are trying to converge.
If your convergence studies have indicated that a large ENCUT is necessary for this system, I would recommend using different potentials, especially the hard h variants, which require a higher plane wave cutoff energy, might be suitable. For all elements, you are using h variants are available.
Cheers, Michael