It depends a bit on the exact kind of calculation that you are running. It's known that some ab initio approaches may have that kind of problem, but I don't know what your software actually does when you do "automatized quantum chemical calculations". A bit more information would be useful.

But anyway, if some computational scheme produces imaginary frequencies (they're not negative, the minus sign is just a stupid misleading flag) for simple and well-behaved systems, then the results are obviously not trustworthy. The simple answer would be: try some other method, and see what sticks. Computational chemistry unfortunately offers only very few hard-and-fast rules about what would be the ideal method for any system if you want to get the calculation done within reasonable time.

Dear Matthias, Thank you for the suggestion. In some cases it did help if I first made a semi-empirical geometry optimization and started the the DFT B3LYP Gaussian 6311 optimization afterwards (which, according to my experience gives more reliable frequencies than HF), followed by force field calculations at the optimized geometry. In some other cases it helped if I switched from Cartesian coordinates to internal coordinates, but always. In some cases even the geometry optimization did not  start at the DFT B3LYP Gaussian 6311 level. Unfortunately the program does not tell too mucha bout the details. Although I started my scientific activity in the 1980'swith semiempirical and ab initio molecular modeling but at that time the programs were much less automated. Now I just kind of play with the program available for me. Now I have to rely on "experimental" computational chemistry.

Well, if your program is interfacing Gaussian, I would recommend to look at the output files that Gaussian creates. Those should of course contain useful info if the program aborts with an error.

Maybe if you can locate the output files, you could upload one with spurious frequencies here so I can have a look and try to reproduce the problems. B3LYP/6-311G (or any derivative of that basis set) should in principle be fine, although there are a few things I would recommend. But I'd like to take a closer look at what exactly you're doing first, if that's not a problem.

Dear Matthias, I'll send you some files privately in a message. Thank you for your help.

Imaginary frequencies are usually a sign of structures that aren't in a fully relaxed ground state, usually for example transition states in reactions paths.

If you get unwanted imaginary frequencies in ground state geometry optimisations, it often can help to tighten the convergence criteria, for example by adding the keywords SCF=tight. (Unless you are already using the most recent version of Gaussian, G16, were this is already the default.)

Dear Simon, thank you for your kind note. It turned out at least in one case that the negative frequency was caused by a wrong assumption: in 4-H-1,2,4-triazol-4-amine I expected an aromatic (planar) structure (similarly to pyrrole), but apparently it is saddle point geometry, not the minimum.

The angular function −ωt θ, has slope −ω, which is named a negative frequency. However, if the function is used as the argument of a cosine operator, the result is not the same from cos(ωt − θ). Thus, sin(−ωt θ) is not the same from sin(ωt − θ π). So any sinusoid can be expressed in names of positive frequencies. The symbol of the underlying phase slope is vague. The vagueness is resolved when the cosine and sine operators can be observed together because cos(ωt-θ) attends sin(ωt- θ) by 1/4 cycle (= π/2 radians) when ω > 0, and lags by 1/4 cycle when ω < 0. Similarly, a vector, (cos t, sin t), rotates counter-clockwise at 1 radian/sec, and completes a circle every 2π seconds, and the vector (cos −t, sin −t) rotates in the other direction.