In this work, we present a gradient-based optimization method to optimize the geometrical properties of metasurfaces based on nano-plasmonic structures, aiming to enhance electric field intensity for applications including high-harmonic generation and surface-enhanced Raman scattering sensing. Our approach involves developing a data-driven deep learning simulator that estimates the electromagnetic response. Specifically, the simulator predicts the electric field distribution at a given cross section of a plasmonic meta-atom based on its geometry. Since the simulator is differentiable, it enables the optimization of various powers of the electric field intensity by coupling it with a generator model and a suitable loss function. With respect to conventional methods that rely on trial and error for optimization and miss atoms intercoupling, our approach systematically considers the collective behavior of the metasurface and efficiently explores the design space. Furthermore, it is capable of capturing multiple local minima within these electrodynamic systems, with optimization results well beyond the training dataset. The generalization capability of the simulator and the behavior of the optimized geometries are validated against a finite element method numerical model.