Rational drug design holds the potential for developing more effective drugs in shorter time frames than traditional trial-and-error approaches. Accordingly, software tools are needed that automate the search process and can computationally direct the changes in ligand geometry and charge distribution to selectively bind at specified receptor sites. Such analysis entails two components: a module that evaluates a cost function (e.g., binding energy/affinity/specificity) for a given parameter set (e.g., charge distribution); and an algorithm for adjusting the parameters to optimize the cost function. Current analysis options include generalized Born (GB), Poisson-Boltzmann (PB) based descriptions and molecular dynamics (MD). GB methods lack geometric detail and MD is expensive so that the balance of modeling fidelity and computational cost in PB methods would seemingly commend them for drug design applications. However, PB methods inadequately characterize the solvent near the moving/vibrating molecular surface, resulting in limited agreement with MD and experimental data; exhibit numerical difficulties (e.g., singular solutions, discontinuous forces); and are overly dependent on choice of surface definition thus limiting predictive reliability. The proposed effort has two core aims. The first is to address the numerical, computational and/or physical modeling deficiencies of PB solvers, improve agreement with MD predictions, and better match experimental data using a statistical description of the moving surface and surrounding ion distributions. An ancillary advantage of this description is that the resulting forces (gradients) are smooth and well-behaved, and thus well-suited for formal optimization-based drug design. The second overall aim pursues this opportunity by formulating and implementing an efficient means of adjusting the design parameters to optimize binding behavior. The approach includes an innovative and computationally efficient means of extracting the required parameter sensitivities from grid-based size-modified PB analysis using an adjoint variable-based method to find all parameter sensitivities in one calculation with computational effort similar to solving the PB Equation. A corollary benefit for both aims is implementation of the algorithms on a novel adaptive Cartesian mesh structure [5, 6] that provides variable mesh spacing, a singularity-free representation of the solution, and fast iterative solution convergence, and has been demonstrated to lower computational resource requirements by over an order of magnitude compared to regular lattice codes. The proposed effort will: (1) Define and incorporate smooth dielectric and ion distributions; (2) benchmark these distributions against MD; (3) validate the approach against well-established small molecule, biomolecular hydration and receptor-ligand binding data; and (4) implement a fast sensitivity analysis using adjoint variables. Successful development will provide a dependable and rational size-modified PB-based molecular design methodology thereby expediting structure-based drug development and significantly improving the commercial viability of these methods.