The electrostatic potential around a molecule is often used to describe reactions, binding, and catalysis mechanisms or to serve as a descriptor in structure-activity relationships and molecular similarity studies. Often, very accurate descriptions of this property are needed that traditionally can be obtained, at least for small molecules, by quantum chemical calculations. The aim of this paper is to extend ab initio-quality quantum chemical accuracy to larger molecules such as proteins. The additive fuzzy density fragmentation (AFDF) principle and the adjustable density matrix assembler (ADMA) method are used to divide large molecules into fuzzy fragments, for which quantum chemical calculations can be done directly using smaller, "custom-made" parent molecules including all the local interactions within a preset distance limit. In the next step, the obtained density matrices of electron density fragments are combined to approximate the global density matrix and the electron density of the whole molecules. These ADMA electron densities are then used to calculate ab inito-quality electrostatic potentials of the large molecules. The accuracy of the method is analyzed in detail by two test cases of a penta- and a hexapetide, and the efficiency of the technique is demonstrated by the calculation of the electrostatic potential of the protein crambin.