When calculating pKa, the main problem is to correctly estimate the solvation energy of all particles involved in the process. You can see a detailed description of how this was done in our work DOI: 10.1007/s00894-018-3726-4 Perhaps this will answer your question.

Please refer to our latest manuscript for the equations;

Article On The Nature of Fe=Oaq in Aqueous Media: A DFT analysis

You should find out the free energy change for the reaction and divide it using 2.303RT.

As mentioned by Alexey N. Masliy accurate usage of solvation energies is vital in the estimation of pKa.

Best wishes.