This is unfortunately not a real answer to your question as such, but I'd like to provide some pointers why this approach is a lot more complicated than it seems at first. It might very well be (currently) impossible.
That is probably also a reason why you don't find a lot of publications on the theory of it and so far no black-box method has been developed. Solvation is still one of the major problems in computational chemistry, hence the use of implicit solvents and counter-charges. Such methods naturally come with the implicit disclaimer that these are approximate corrections only.
The most elaborate method I know so far is COSMO-RS, which still has large error margins, at least that is the my current state of knowledge.
I think for calculations you would at least need to run molecular dynamics simulations with many explicit solvent molecules to describe the equilibrium. This is important, because the proton which dissociates will be caught in the bulk of water molecules and a highly complex network of hydrogen bonds. This is very well demonstrated by the Grotthuss mechanism, which shows how easily bond lengths can change. I would say it is impossible to describe this complex system with a single quantum chemical calculation.
Furthermore, DFT is notoriously bad at describing non-covalent interactions, and including dispersion corrections will only get you so far. I have not checked recently, but I believe there are quite a few publications discussing various sizes of water clusters on this level. I have neither used nor heard of AFPD before, so it would be wise to thoroughly search for literature (benchmarks) about its performance on these clusters. Suffice to say, that one particular functional will not give you a reliable result as you have nothing to compare it to. Especially in such complicated cases you have to calibrate your results.
Another problem here is that you are the describing heterolytic bond dissociation. More often than not this may lead to additional problems. I don't always have confidence in my results when I treat simple donor-acceptor complexes where there is no charge separation and all singlet species. A supermolecular approach would probably a better choice, but could heavily suffer from a BSSE. However, it might be worth considering $$\ce{H2O\bond{...}HA <=> H2OH\bond{...}A}.$$
I think the key for calculating acid dissociation constants is treating the whole system as one ensemble.
Recall the common definition of the equilibrium constant for the reaction
\begin{align}\ce{AH (aq) &<=> A- (aq) + H+ (aq)},&
K &= \frac{a(\ce{A-})\,a(\ce{H+})}{a(\ce{HA})}.\end{align}
We conveniently ignore water here, because it is abundant and we assume ideal solution. However, when your system is too small, you will get a large error for $K$. This is most likely due to your solvated proton, which will be described worst with solvation models. Explicitly you are writing
$$\ce{HA (aq) + H2O (aq) -> A− (aq) + H3O+ (aq)},$$
which is not really what you are calculating. Instead you are (more ore less) using already approximations to treat thermal and solvent corrections
$$\ce{HA (g, th, solv) + H2O (g, th, solv) -> A− (g, th, solv) + H3O+ (g, th, solv)}.$$
And this doesn't even consider that you are calculating two states in an equilibrium and you will get the equilibrium constant (using activities, not concentrations) for only these two states (which might not even be the ground sates),
$$K = \exp\left\{-\frac{\Delta G}{RT}\right\}.$$
You have also not accounted for the possible conformational space of $\ce{HA}$/$\ce{A-}$, which might become very important as soon as you have ethyl moieties, or even cyclohexyl moieties.
These are just some issues that immediately come to mind, there might very well be more underlying insufficiencies. I have looked into this a couple of years ago and came to the conclusion not to bother with it any more.
What can you do instead?
Assessing acid strengths does not need to be at the level of what experimental chemists understand. If they really need the $\ce{pKa}$ they should run the appropriate measurements and obtain reliable results instead of educated guessing. If they want you to classify whether an acid is stronger or weaker, use a rigorous approach to obtain the conformational space for all species and compare them with isodesmic equations. This will have the advantage that insufficiencies of the methodology will likely cancel each other.
Another purely theoretical approach would be calculating gas-phase proton affinities. From this you can judge if a compound behaves stronger or weaker. Keep in mind that this approach is better suited for bases or nucleophiles. If you have the resources, you can also evaluate Lewis acidities,
$$\ce{Base + BX3 <=> Base-BX3}.$$
This should give you at least a good idea how on a relative level your compounds behave.
You can also use a first-principle derived value for your proton, which will likely cause the largest error. A quick check turned up a few promising results. (For example: Chang-Guo Zhan and David A. Dixon. J. Phys. Chem. A 2001, 105 (51), 11534–11540.) You can then go ahead and calculate know compounds and try to find a correlation between $\ce{pKa}$ and $\Delta G$. This way of "calibration" might lead to some significant error compensation.
I would also recommend to look at COSMO-RS in more detail. There might be a lot of development that I have missed and they might have found a way to obtain more reliable results. As far as I remember, the foundation of this was solid and a good step into easier capturing solvent effects without the need to resort to full blown MD simulations.
