The evaluation of molecular electrostatic potential (ESP) is a performance bottleneck for many computational chemical tasks like restrained ESP charge fitting or quantum mechanics/molecular mechanics simulations. In this paper, an efficient algorithm for the evaluation of ESP is proposed. It regroups the expression in terms of primitive Gaussian type orbitals (GTOs) with identical angular momentum types and nuclei centers. Each term is calculated using a computerized optimized code. This algorithm was integrated into the wavefunction analysis program Multiwfn and was tested on several large systems. In the cases of dopamine and remdesivir, the performance of this algorithm was comparable to or better than some popular state-of-the-art codes. For meta1-organic framework-5, where the number of GTOs and ESP points is 4840 and 259 262, respectively, our code could finish the evaluation in 1874 seconds on ordinary hardware. It also exhibits good parallelization scaling. The source code of this algorithm is freely available and can become a useful tool for computational chemists. The computerized optimized code enables Multiwfn to carry out large scale electronstatic potential analyses for systems of considerable size.