  A computational procedure is developed for the efficient calculation of
derivatives of integrals over non-separable Gaussian-type basis functions, used
for the evaluation of gradients of the total energy in quantum-mechanical
simulations. The approach, based on symbolic computation with computer algebra
systems and automated generation of optimized subroutines, takes full advantage
of sparsity and is here applied to first energy derivatives with respect to
nuclear displacements and lattice parameters of molecules and materials. The
implementation in the \textsc{Crystal} code is presented and the considerably
improved computational efficiency over the previous implementation is
illustrated. To this purpose, three different tasks involving the use of
analytical forces are considered: i) geometry optimization; ii) harmonic
frequency calculation; iii) elastic tensor calculation. Three test case
materials are selected as representatives of different classes: i) a metallic
2D model of the Cu (111) surface; ii) a wide-gap semiconductor ZnO crystal,
with a wurtzite-type structure; and iii) a porous metal-organic crystal, namely
the ZIF-8 Zinc-imidazolate framework. Finally, it is argued that the present
symbolic approach is particularly amenable to generalizations, and its
potential application to other derivatives is sketched.
