CROSS-REFERENCES TO RELATED APPLICATIONS The present application claims priority from and is a nonprovisional application of J.S. Provisional Application No. 60/511,477, entitled "Method and Apparatus for Assessing I Likelihood of Combination of Molecules Using a Basis Expansion" filed October 14,2003, he entire contents of which are herein incorporated by reference for all purposes. FIELD OF THE INVENTION The present invention generally relates to bioinformatics, proteomics, molecular nodeling, computer-aided molecular design (CAMD), and more specifically computer-aided Irug design (CADD) and computational modeling of molecular combinations. BACKGROUND OF THE INVENTION An explanation of conventional drug discovery processes and their limitations is useful for understanding the present invention. Discovering a new drug to treat or cure some biological condition, is a lengthy and expulsive process, typically taking on average 12 years and $800 million per drug, and taking possibly up to 15 years or more and $ 1 billion to complete in some cases A goal of a drug discovery process is to identify and characterize a chemical compound or ligand biomolecule, that affects the function of one or more other biomolecules (i.e., a drug "target") in an organism, usually a biopolymer, via a potential molecular interaction or combination. Herein fee term biopolymer refers to a macromolecule that comprises one or more of a protein, nucleic acid (DNA or RNA), peptide or nucleotide sequence or any portions or fragments thereof Herein the term biomolecule refers to a chemical entity that comprises one or more of a biopolymer, carbohydrate, hormone, or other molecule or chemical compound, either inorganic or organic, including, but not limited to, synthetic, medicinal, drug-like, or natural compounds, or any portions or fragments thereof The target molecule is typically a disease-related target protein or nucleic acid for which it is desired to affect a change in function, structure, and/or chemical activity in order id in the treatment of a patient disease or other disorder. In other cases, the target is a molecule found in a disease-causing organism, such as a virus, bacteria, or parasite, that en affected by the drug will affect the survival or activity of the infectious organism. In other cases, the target is a biomolecule of a defective or harmful cell such as a cancer cell, yet other cases the target is an antigen or other environmental chemical agent that may luce an allergic reaction or other undesired immunological or biological response. The ligand is typically a small molecule drug or chemical compound with desired iig-like properties in terms of potency, low toxicity, membrane permeability, solubility, lemical / metabolic stability, etc. In other cases, the ligand may be biologic such as an jected proteirrbased or peptide-based drug or even another full-fledged protein. In yet iiM-cases the ligand may be a chemical substrate of a target enzyme. the ligand may even e covalently bound to the target or may in fact be a portion of Ae protein, e.g., protein scondary structure component, protein domain containing or near an active site, protein ubunit of an appropriate protein quaternary structure, etc. Throughout the remainder of the background discussion, unless otherwise specifically lifferentiated, a (potential) molecular combination will feature one ligand and one target, the igand and target will be separate chemical entities, and the ligand will be assumed to be a sheroical compound while the target will be a biological protein (mutant or wild type). Note that the frequency of nucleic acids (both DNA/RNA) as targets will likely increase in coming years as advances in gene therapy and pathogenic microbiology progress. Also the term "molecular complex" will refer to the bound state between the target and ligand when interacting with one another in the midst of a suitable (often aqueous) environment. A "potential" molecular complex refers to a bound state that may occur albeit with low probability and therefore may or may not actually form under normal conditions. The drug discovery process itself typically includes four different subprocesses: (1) target validation; (2) lead generation / optimization; (3) preclinical testing; and (4) clinical trials and approval. Target validation includes determination of one or more targets that have disease relevance and usually takes two-and-a-half years to complete. Results of the target validation phase might include a determination that the presence or action of the target molecule in an organism causes or influences some effect that initiates, exacerbates, or contributes to a disease for which a cure or treatment is sought. In some cases a natural binder or substrate for the target may also be determined via experimental methods. Lead generation typically involves the identification of lead compounds that can bind ) the target molecule and thereby alter the effects of the target through either activation, eactivation, catalysis, or inhibition of the function of the target, in which case the lead would e a viewed as a suitable candidate ligand to be used in the drag application process. Lead •ptimization involves the chemical and structural refinement of lead candidates into drag trecursors in order to improve binding aflanity to the desired target, increase selectivity, and ilso to address basic issues of toxicity, solubility, and metabolism. Together lead generation and lead optimization typically takes about three years to complete and might result in one or nore chemically distinct leads for further consideration. In preclinical testing, biochemical assays and animal models are used to test the selected leads for various pharmacokinetic factors related to drug absorption, distribution, metabolism, excretion, toxicity, side effects, and required dosages. This preclinical testing takes approximately one year. After the preclinical testing period, clinical trials and approval take another six to eight or more years during which the drug candidates are tested on htmian subjects for safety and efficacy. Rational drug design generally uses structural information about drug targets (structure-based) and/or their natural ligands (ligand-based) as a basis for the design of effective lead candidate generation and optimization. Structure-based rational drug design generally utilizes a three-dimensional model of the structure for the target. For target proteins or nucleic acids such structures may be as the result of X-ray crystallography / NMR or other measurement procedures or may result from homology modeling, analysis of protein motifs and conserved domains, and/or computational modeling of protein folding or the nucleic acid equivalent. Model-built structures are often all that is available when considering many membrane-associated target proteins, e.g., GPCRs and iorrcharmels. The structure of a ligand may be generated in a similar marmer or may instead be constructed ab initio from a known 2-D chemical representation using fundamental physics and chemistry principles, provided the ligand is not a biopolymer. Rational drug design may incorporate the use of any of a number of computational components ranging from computational modeling of target-ligand molecular interactions and combinations to lead optimization to computational prediction of desired drug-like biological properties. The use of computational modeling in the context of rational drug design has been largely motivated by a desire to both reduce the required time and to improve the focus and efficiency of drug research and development, by avoiding often time consuming and costly efforts in biological "wet" lab testing and the like. Computational modeling of target-ligand molecular combinations in the context of ad generation may involve the large-scale irrsilico screening of compound libraries (i.e., jrary screening), whether the libraries are virtually generated and stored as one or more >mpound structural databases or constructed via combinatorial chemistry and organic /nthesis, using computational methods to rank a selected subset of ligands based on imputational prediction of bioactivity (or an equivalent measure) with respect to the itended target molecule. Throughout the text, the term "binding mode" refers to the 3-D molecular structure of potential molecular complex in a bound state at or near a minimum of the binding energy i.e., maximum of the binding affinity), where the term "binding Oiergy" (sometimes aterchanged with "binding affinity*' or "binding free energy") refers to the change in free snergy of a molecular system upon formation of a potential molecular complex, i.e., the ransition from an unbound to a (potential) bound state for the ligand and target. Here the erm free energy generally refers to both enthalpic and entropic effects as the result of )hysical interactions between the constituent atoms and bonds of the molecules between hemselves (i.e., both intermolecular and intramolecular interactions) and with their surrounding environment. Examples of the free energy are the Oibbs free energy encountered in the canonical or grand canonical ensembles of equilibrium statistical mechanics. In general, the optimal binding free energy of a given target-ligand pair directly correlates to the likelihood of formation of a potential molecular complex between the two molecules in chemical equilibrium, though, in truth, the binding free energy describes an ensemble of (putative) complexed structures and not one single binding mode. However, in computational modeling it is usually assumed that the change in free energy is dominated by a single structure corresponding to a minimal energy. This is certainly true for tight binders (pK ~ 0.1 to 10 nanomolar) but questionable for weak ones (pK ~ 10 to 100 micromolar). The dominating structure is usually taken to be the binding mode. In some cases it may be necessary to consider more than one alternative-binding mode when the associated system states are nearly degenerate in terms of energy. It is desirable in the drug discovery process to identify quickly and efficiently the optimal dockmg configurations, i.e., bindmg modes, of two molecules or parts of molecules. Efficiency is especially relevant in the lead generation and lead optimization stages for a drug discovery pipeline, where it is desirable to accurately predict the binding mode for possibly illions of potential molecular complexes, before submitting promising candidates to further lalysis. Binding modes and binding affinity are of direct interest to drug discovery and itional drug design because they often help indicate how wdl a potential drug candidate may jrve its purpose. Furthermore, where the binding mode is determinable, the action of the rug on the target can be better understood. Such understanding may be useful when, for sample, it is desirable to further modify one or more characteristics of the ligand so as to nprove its potency (with respect to the target), binding specificity (with respect to other irgets), or other chemical and metabolic properties. A number of laboratory methods exist for measuring or estimating affinity between a arget molecule and a ligand. Often the target might be first isolated and then mixed with the igand in vitro and the molecular interaction assessed experimentally such as in the myriad nochemical and functional assays associated with higjh throughput screening. However, such neihods are most useful where the target is simple to isolate, the ligand is simple to nanufacture and the molecular interaction easily measured, but is more problematic when the target carmot be easily isolated, isolation interferes with the biological process or disease patheray, the ligand is difficult to synthesize in sufficient quantity, or where the particular target or ligand is not well characterized ahead of time. In the latter case, many thousands or millions of experiments might be needed for all possible combinations of the target and ligands, making the use of laboratory methods unfeasible. While a number of attempts have been made to resolve this bottleneck by first using specialized knowledge of various chemical and biological properties of the target (or even related targets such as protein family members) and/or one or more already known natural binders or substrates to the target, to reduce the number of combinations required for lab processing, this is still impractical and too expensive in most cases. Instead of actually combining molecules in a laboratory setting and measuring experimental results, another approach is to use computers to simulate or characterize molecular interactions between two or more molecules (i.e., molecular combinations modeled in silico). The use of computational methods to assess molecular combinations and interactions is usually associated with one or more stages of rational drug design, whether structure-based, ligand-based, or both. The computational prediction of one or more binding modes and/or the computational assessment of the nature of a molecular combination and the likelihood of formation of a potential molecular complex is generally associated with the term "docking" in the art To e, conventional "docking" methods have included a wide variety of computational hniques as described in the forthcoming section entitled "REFERENCES & PRIOR ART". Whatever the choice of computational docking method there are inherent trade-offs tween the computational complexity of both the underlying molecular models and the linsic numerical algorithms, and the amount of computing resources (time, number of *Us, number of simulations) that must be allocated to process each molecular combination, >r example, while highly sophisticated molecular dynanaics simulations (MD) of the two olecules surrounded by explicit water molecules and evolved over trillions of time steps lay lead to higher accuracy in modeling the potential molecular combination, the resultant smputational cost (i.e., time and computing power) is so enormous that such simulations are vtractable for use with more than just a few molecular combinations. One major distinction amongst docking methods as applied to computational nodeling of molecular combinations is whether the ligand and target structures remain rigid hroughout the course of the simulation (i.e., rigid-body docking) vs. the ligand and/or target >eing allowed to change their molecular conformations (i.e., flexible docking). In general, :he latter scenario involves more computational complexity, though flexible docking may often achieve higher accuracy than rigid-body docking when modeling various molecular combinations. That being said rigid-body docking can provide valuable insight into the nature of a molecular combination and/or the likelihood of formation of apotential molecular complex and has many potential uses within the context of rational drug discovery. For instance rigid-body docking may be appropriate for docking small, rigid molecules (or molecular fragments) to a simple protein with a well-defined, nearly rigid active site. As another example, ri^d-body docking may also be used to more efficaently and rapidly screen out a subset of likely nonactive ligands in a molecule library for a given target, and then applying more onerous flexible docking procedures to the surviving candidate molecules. Rigid-body docking may also be suitable for de novo ligand design and combinatorial library design. Moreover, in order to better predict the binding mode and better assess the nature and/or likelihood of a molecular combination when one or both molecules are likely to be flexible, rigid-body docking can be used in conjunction with a process for gaierating likely yet distinct molecular conformers of one or both molecules for straightforward and efBcient virtual screening of a molecule library against a target molecule. However, as will be discussed, even rigid body docking of molecular combinations can be computationally expensive and thus there is a clear need for better and more efficient computational methods ed on rigid body docking when assessing the nature and/or likelihood of molecular nbinations. As outlined in the section entitled "REFERENCES & PRIOR ART', conventional mputational methods for predicting binding modes and assessing the nature and/or :elihood of molecular combinations in the context of rigid-body docking include a wide riety of techniques. These include methods based on pattern matcdiing (often graph-based), aximization of shape complementarity (i.e., shape correlations), geometric hashing, pose ustering, and even the use of one or more flexible docking methods with the simplifying m-time condition that both molecules are rigid. Of special interest to this invention is class of rigid-body doddng techniques based on he maximization of shape complementarity via evaluation of the spatial correlation between wo representative molecular surfaces at different relative positions and orientations. One sxample is the "Hex" docking software described in Ritchie, D. W. and Kemp. G. J. L, 'Protein Docking Using Spherical Polar Fourier Correlations", (2000), Proteins: Structure. Function, and Genetics, 39,178-194; (hereinafter, "Ritchie et ar"), all of which is hereby incorporated by reference in their entirety. Further examples include the "FTDOCK" docking software of the Cambridge Crystallographic Data Center described in Aloy, P., Moont, 0., Gabb, H. A., Querol, £., Aviles, F. X., and Sternberg, M. J. E., "Modeling Protein Docking using Shape Complementarity, Electrostatics and Biochemical Information," (1998), Proteins: Structure, Function, and Genetics, 33(4) 535-549; all of which is hereby incorporated by reference in their entirety. Such shape complementarity based methods while typically treating molecules as rigid and thus perhaps less rigorous than their flexible docking counterparts, especially in the context of flexible molecules, is still potentially valuable for the fast, efficient screening of two molecules in ordw to make a preliminary assessment of the nature and/or likelihood of formation of a potential molecular complex of the two molecailes or to make an initial prediction of the preferred binding mode for the molecular combination. Such a preliminary assessment may significantly reduce the number of candidates that must be further screened in silico by another more computationally costly docking method. Moreover, the utility of computing shape complementarity has been demonstrated with respect to multiple proteirr protdn systems, including both enzyme-inhibitor and antibody-antigen, as per FTDOCK and 'Ritchie et al. Previous formulations for the computation of shape complementarity generally fall tto four categories. The first category involves the calculation of a spatial correlation in the spatial omain between two volumetric functions describing a representative molecular surface for ach molecule, where the spatial correlation between two 3-D complex functions, f(r) and l(f), is calculated as follows: f vhere f denotes the complex conjugate and * denotes convolution. In the spatial domain, the spatial correlation is performed by converting the integrals into summations and directly jomputing over a sampling space comprising three translation variables with a specified sampling resolution. A search is then conducted by reevaluating the spatial correlation for sach new and different relative orientation of the molecules. Those configurations that show the highest net spatial correlation are typically selected as possible candidate binding modes. However, the method of directly computing the spatial correlation in the spatial domain is often very computationally intensive, since if the sampling translation grid is a M x M X M grid, the above spatial correlation calculation requires 0(M*) operations. For instance, when M= 256, there are more than 2.8 x lO'* multiplication operations required. Furthermore, the 0(M*) calculations must be performed for every relative orientation of the two molecules, making the total number of calculations impractical at best. The second category involves the calculation of a spatial correlation in the frequency domam between two volumetric functionss describing a representative molecular surface for each molecule, where the spatial correlation betweeen two 3-D functions, f and g, is still defined as before, but a firequency space decomposition, such as a Fourier transform, is used in order to reduce the number of calculations. For a full description of the Fourier transform refer to Press, W. H., Flarmery, B. P., Teukolsky, S. A., and Vetteriing, W. T., "Numerical Recipes in C: The Art of Scientific Computing", Cambridge University Press (1993) (hereinafter, "Press et af'), hereby incorporated by reference in its entirety. To compute a 3-D spatial correlation in Fourier space, one can use the following relation, also known as the Correlation Theorem (Press et at). The convolution of two complex 3-D functions f(?) with g(r) is given by: f*g = F-'(F(u)G{0)) [Eqn.2] ere F' refers to the inverse Discrete Fourier Transform (Press et at), F(U) is Discrete irier Transform of the complex conjugate of f( r) and G(U) is the Discrete Fourier insform of g( 7). Similar formulae can be generated for other frequency decompositions sides Fourier, such as a Laplace transform, a Disca-ete Cosine transform, and others, as scribed in Arfken, G. B., and Weber H. J., "Mathematical Method for Physicist", Harcourt tcademic Press (2000), (hereinafter, "Arfiken et aP'), hereby incorporated by reference in its tirety. For the same M x M x M grid, the frequency based evaluation of the spatial •rrelation will require approximately 0(3M^ ln(M)) operations where ln(M) denotes the iturallogaridimofM. While the number ofoperations decreases substantially wh«i the FTs are used as opposed to direct computation in the spatial domain, the amount of memory orage and/or the amount of data that must be read from storage must still be taken into ;count, i.e., the input/output (I/O) bandwidth requirement. For example, for M = 256 at 16 bit precision, 800 Mbits are required for computing le 3-D spatial correlation using DFTs for just one relative orientation. Generally, this is a ery large amount of data for storage directly in memory and would require millions of clock ycles to fetch from one or more DRAM chips with current DRAM and I/O bus technology, vloreover, the 0(3M^ ln(M)) calculations and the access of hundreds of millions of data bits nust be performed for every relative orientation ofthe two molecules, making the entire process onerous when considering possibly millions of such relative orientations in the course Df a high resolution search of fee shape complementarity space. For this reason, Fourier based methods for evaluating shape complementarity often take hours on conventional compute software in order to complete for large protein systems, for instance, as in FTDOCK, and as such are not suitable for large-scale screening. The third category involves the least-square minimization (or equivalent minimization) of separation distances between critical surface and/or fitting points that represent the molecular surfaces of the two molecules. Examples include the Patchdock docking software written by Nussinov-Wolfson Structural Bioinformatics Group at Tel-A\'iv University, based on principles described in Lin, S. L., Nussinov, R., Fischer, D,, and Wolfson, H. J., "Molecular surface representations by sparse critical points", (1994) Proteins: Structure, Function, and Genetics, 18,94-101; all of which are hereby incorporated by reference in their entirety. Such methods often suffer from degraded accuracy, especially when the molecular jrface geometry is complex or when the ligand molecule is very small relative to the protein sceptor and/or characterized by poor binding affinities. Moreover, the cost of computing the arface critical points is often itself very expensive. The number of computations associated with the three method categories described bove renders the process impractical for use with conventional computer software and ardware configurations when performing large-scale screening. Moreover, the above lethods are not practical for high accuracy prediction of the binding mode due to the equirement of a high resolution of the associated sampling space, A fourth category has been developed for the efficient estimation of shape iomplementarity based on the decomposition of two volumetric functions describing a epresentative molecular surface for each molecule onto an appropriate orthogonal basis set, luch as a radial-spherical harmonics expansion, as described in Ritchie et al. The chief idvantage of this type of method is that the required number of calculations scale linearly . with the desired number of sampled configurations, thus allowing for a dense sampling of the geometric shape complementarity. Moreover, the computing time is roughly invariant vn&i respect to the sizes of the two molecules and is thus suitable for proteirrprotein doddng. However, to achieve high accuracy for complex molecular surface geometries, it is necessary to perform the orthogonal basis expansion with a large expansion order and as such the total computing time can be quite large. Furthermore, current methods such as those outlined in Ritchie et al are not amenable to implementation in customized or other application specific hardware for use in large-scale screening. In summary, it is desirable in the drug discovery process to identify quickly and efficiently the optimal configurations, i.e., binding modes, of two molecules or parts of molecules. Efficiency is especially relevant in the lead generation and lead optimization stages for a drug discovery pipeline, where it may be desirable to accurately predict the binding mode and binding affinity for possibly millions of potential target-ligand molecular combinations, before submitting promising candidates to further analysis. There is a clear need then to have more efficient systems and methods for computational modeling of the molecular combinations with reasonable accuracy. In general, the present invention relates to an efficient computational method for analysis of molecular combinations based on maximization of shape complementarity over a set of configurations of a molecular combination through computation of a basis expansion representing molecular shapes in a coordinate system. Here the analysis of the molecular nbination may involve the prediction of likelihood of formation of a potential molecular uplex, the prediction of the binding mode (or even additional alternative modes) for the mbination, the characterization of the nature of the interaction or binding of various mponents of the molecular combination, or even an approximation of binding affinity for s molecular combination based on a shape complementarity score or an equivalent sasure. The invention also addresses and solves the various hardware implementation irdles and bottlenecks associated with current conventional methods. REFERENCES & PRIOR ART Prior art in the field of the current invention is heavily documented: the following ies to summarize it. . Drews [1] provides a good overview of the current state of drug discovery. In [2] bagyan and Totrov show the state of high throughput docking and scoring and its pplications. Lamb et al [3] further teach a general approach to the design, docking, and irtual screening of multiple combinatorial libraries against a family of proteins, finally Vaskowycz et al [4] describe the use of multiple computers to accelerate virtual screening of large ligand library against a specific target by assigning groups of ligands to specific omputers. [1] J. Drews, "Drug Discovery: A Historical perspective," Science 287,1960- 1964(2000). [2] Ruben Abagyan and Maxim Totrov, "High-throughput docking for lead generation". Current Opinion in Oiemical Biology 2001,5:375-382. [3] Lamb, M. L.; Burdick, K. W.; Toba, S.; Young, M. M.; Skillman, A. G. et al "Design, docking, and evaluation of multiple libraries against multiple targets". Protdns 2001,42,296-318. [4] Waszkowycz, B., Perkins, T.D.J., Sykes, R.A., Li, J. "Large-scale virtual screening for discovering leads in the post-gaiomic «a", IBM Systems Joumal, Vol. 40, No. 2 (2001). There are a number of examples of software tools currently used to perform docking simulations. These methods involve a wide range of computational techniques, including use of a) rigid-body patterrrmatching algorithms, either based on surface correlations, use of geometric hashing, pose clustering, or graph patterrrmatching; b) fragmental-based methods, including incremental construction or "place and join" operators; c) stochastic optimization thods including use of Monte Carlo, simulated armealing, or genetic (or memetic) prithms; d) molecular dynamics simulations; or e) hybrids strategies derived thereof. The earliest docking software tool was a graph-based rigid-body patterrrmatching jorithm called DOCK [5,6, 7], developed at UCSF back in 1982 (vl .0) and now up to v5.0 ifh extensions to include incremental construction). Other examples of graph-based tterrrmatching algorithms include CLDC [8] (which in turn uses GRID [9]), FLOG [ 10] and GIN [11]. [5] Shoichet, B.K., Bodian, D.L. and Kuntz, ID., "Molecular docking using shape descriptors",J.Co7Mp.C/jgm., Vol. 13 No. 3,380-397 (1992). , [6] Meng, E.G., Gschwend, D.A., Blaney, J.M., and I.D. Kuntz, "Orientational sampling and rigid-body minimization in molecular docking", Proteins: Structure, Function, and Genetics. Vol. 17,266-278 (1993). [7] Ewing, T. J. A. and Kuntz, I. D., "Critical Evaluation of Search Algorithms for Automated Molecular Docking and Database Screening", J. Computational aiemisrrv. Vol! 8 No. 9,1175-1189 (1997). [8] Lawrmce, M.C. and Davis, P.C.; "CLDC: A Search Algorithm for Finding Novel Ligands Capable of Binding Proteins of Known Three-Dimensional Structure", Proteins, Vol. 12,31-41 (1992). [9] Kastenholz, M. A., Pastor, M., Crudani, G., Haaksma, E. E. J., Fox, T., "GRID/CPCA: A new computational tool to design selective ligands", J. Medicinal Chemistry, Vol. 43, 3033-3044 (2000). [10] Miller, M. D., Kearsley, S. K., Underwood, D. J, and Sheridan, R. P., "FLOG: a system to select *qasrflexible'ligands complementary to a receptor of known tliree-dimensional structure",/ Computer-Aided Molecular Design^ Vol. 8 No.2,153-174 (1994). [11] Sobolev, v., Wade, R. C, Vriend, G. and Edelman, M., "Molecular docking using surface complementarity", Proteins, Vol. 25,120-129 (1996). Other rigid-body patterrrmatching docking software tools include the shape-based correlation methods of FTDOCK [12] and HEX [13], the geometric hashing of Fischer et al [14], or the pose clustering of Rarey et al [15]. [12] Katchalsld-Katzir, E., Shariv, I., Eisenstein, M., Friesem, A. A., Aflalo, C, and Vakser, I.A., "Molecular siirface recognition: Determination of geometric fit between proteins and their ligands by correlation techniques", Proceedings of the Natiojial Acadeyny of Sciences of the United States of America, Vol. 89 No. 6,2195-2199 (1992). [13] Ritchie, D. W. and Kemp. G. J. L., "Fast Computation, Rotation, and Comparison of Low Resolution Spherical Harmonic Molecular Surfaces", / Computational Chemistry, Vol 20 "No. 4,3S3-395 (1999). [ 14] Fischer, D., Norel, R., Wolfson, H. and Nussinov, R., "Surface motifs by a computer vision technique: searches, detection, and implications for proteirr ligand recognition", Proreirw, Vol. 16, 278-292 (1993). [15] Rarey,M.,Wefing,S., and Lengauer,T., "Placement of mwiium-sized molecular fragments into active sites of proteins", J. Computer-Aided Molecular Design, Vol. 10,41-54 (1996). In general, rigid-body patterrrmatching algorithms assume that both the target and igand are rigid (i.e., not flexible) and hence may be appropriate for docking small, rigid nolecules (or molecular fragments) to a simple protein with a well-defined, nearly rigid ictivesite. Thus this class ofdocking tools may be suitable for de novo ligand design, combinatorial library design, or straightforward rigid-body screening of a molecule library containing multiple conformers per ligand. Incremental construction based docking software tools include FlexX [16,17] from Tripos (licensed from EMBL), Hammerhead [18], DOCK v4.0 [7] (as an option), and the nongreedy, backtracking algorithm of Leach et al [19]. Programs using incremental construction in the context of cfe novo ligand design include LUDI [20] (from Accelrys) and GrowMol [21]. Docking software tools based on 'place and join' strategies include DesJarlaise/d!/[22]. [16] Kramer, B., Rarey, M. and Lengauer, T., "Evaluation of the FlexX incremental construction algorithm for proteirrligand docking", Proteins, Vol. 37,228-241(1999). [17] Rarey, M., Kramer, B., Lengauer, T., and Klebe, G.," A Fast Flexible Docking Method Using An Incremental Construction Algorithm", J. Mol. 5io/., Vol. 261,470-489 (1996). [18] Welch, W., Ruppert, J. and Jain, A. N., "Hammerhead: Fast, fully automated docking of flexible ligands to protein binding sites", Chemical Biology, Vol. 3, 449-462(1996). [19] Leach, A.R., Kuntz, I.D., "Conformational Analysis of Flexible Ligands in Macromolecular Receptor Sites", J. Comp. Chem., Vol. 13,730-748 (1992). [20] Bohm,H. J., "the computer program LUDI: a new method for the de novo design of enzyme inhibitors", J. Computer-Aided Molecular Design, Vol. 6, 61-78(1992). [21] Bohacek, R. S. and McMartin, C, "Multiple Hig^y Diverse Structures Complementary to Enzyme Binding Sites: Results of Extensive Application of a de Novo Design Method Incorporating Combinatorial Growth", J. American Chemical Society,Vol 116, 5560-5571 (1994). [22] DesJarlais, R,L., Sheridan, R.P., Dixon, J.S., Kuntz, LD., and Venkataraghavan, R., "Docking Flexible Ligands to Macromolecular Receptors by Molecular Shape", J. Med C/iem., Vol. 29,2149-2153 (1986). Incremental construction algorithms may be used to model docking of flexible ligands 0 a rigid target molecule with a well-charactmzed active site. They may be used wh«i screening a library of flexible ligands against one or more targets. They are often x)mparatively less compute-intensive, yet consequently less accurate, than many of their stochastic optimization based competitors. However, even FlexX may take on order of < 1-2 minutes to process one target-ligand combination and thus may still be computationally onerous depending on the size of the library (e.g., tens of millions or more compounds). Recently FlexX was extended to FlexE [23] to attempt to account for partial flexibility of the target molecule's active site via use of user-defined ensembles of certain active site rotamers. [23] Claussen, H., Buning, C, Rarey, M., and Lengauer, T., "FlexE: Efficient Molecular Docking Considering Protein Structure Variations", J. Molecular ^/o/og)^ Vol. 308,377-395 (2001). Computational docking software tools based on stochastic optimization include ICM [24] (from MoISoft), GLIDE [25] (from Schrodinger), and LigandFit [26] (from Accelrys), all based on modified Monte Carlo techniques," and AutoDock v.2.5 [27] (from Scripps Institute) based on simulated armealing. Others based on genetic or memetic algorithms include GOLD [28,29], DARWIN [30], and AutoDock v.3.0 [31] (also from Scripps). [24] Abagyan, R.A., Totrov, M.M., and Kuznetsov, D.N., "Biased probability Monte Carlo conformational searches and electrostatic calculations for peptides and proteins", J. Comp. Chem., Vol. 15,488-506 (1994). [25] Halgren. T.A., Murphy, R.B., Friesner, R.A., Beard, H.S., Frye, L.L., Pollard, W.T., and Banks, J.L., "Glide: a new approach for rapid, accurate docking and scoring. 2. Enrichment factors in database screening", JMerfCAewi., Vol. 47 No. 7,1750-1759, (2004). [26] Luty, B. A., Wasserman, Z. R., Stouten, P. F. W., Hodge, C. N., Zacharias, M., and McCammon, J. A., "Molecular Mechanic^Qrid Method for the Evaluation of Ligand-Receptor Interactions", J, Comp. Chem., Vol.16,454-464(1995). [27] Goodsell,D. S. and Olson, A. J., "Automated Docking of Substrates to Proteins by Simulated Aimealing", Proteins: Structure, Function, and Ge7je/ic5, Vol. 8,195-202 (1990). [28] Jones, G., Willett, P. and Glen, R. C, "Molecular Recognition of Receptor Sites using a Genetic Algorithm with a Description of Desolvation", J. Mol. J?w/., Vol. 245,43-53 (1995). [29] Jones, G., Willett, P., Glen, R. C, Leach, A., and Taylor, R., "Development and Validation of a Genetic Algorithm for Flexible Docking", J. Mol Biol., Vol. 267,727-748 (1997). [30] Taylor, J.S. and Burnett, R. M., Proteins, Vol. 41,173-191 (2000). [31] Morris, G. M., Goodsell, D. S., Halliday, R. S., Huey, R., Hart, W. E., Belew, R. K. and Olson, A. J., "Automated Docking Using a Lamarckian Genetic Algorithm and an Empirical Binding Free Energy Function", J. Comp. Chem., Vol. 19,1639-1662(1998). Stochastic optimizatiorrbased methods may be used to model docking of flexible ligands to a target molecule. They generally use a molecular-mechanics based formulation of the affinity function and employ various strategies to search for one or more favorable system energy minima. They are often more compute-intensive, yet also more robust, than their incremental construction competitors. As they are stochastic in nature, different runs or simulations may often result in different predictions. Traditionally most docking software tools using stochastic optimization assume the target to be nearly rigid (i.e., hydrogen bond donor and acceptor groups in the active site may rotate), since otherwise the combmatorial jmplexity increases rapidly making the problem difficult to robustly solve in reasonable me. Molecular dynamics simulations have also been used in the context of computational lodeling of target-ligand combinations. This includes the implementations presented in Di lola et al[32] and Luty et al[16] (along with Monte Carlo). In principle, molecular ynamics simulations inay be able to model protein flexibility to an arbitrary degree. On the ther hand, they may also require evaluation of many fine-grained, time steps and are thus iften very time-consuming (one order of hours or even days pw target-ligand combination). Tiey also often require user-interaction for selection of valid trajectories. Use of molecular lynamics simulations in lead discovery is therefore more suited to local minimization of >redicted complexes featuring a small number of promising lead candidates. [32] Di Nola, A., Berendsen, H. J. C, and Roccatano, D., "Molecular Dynamics Simulation of the Docking of Substrates to Proteins", Proteins, Vol. 19,174-182(1994). Hybrid methods may involve use of rigid-body pattern matching techniques for fast screening of selected low-energy ligand conformations, followed by Monte Carlo torsional optimization of surviving poses, and finally even molecular dynamics refinement of a few choice ligand structures in combination with a (potentially) flexible protein active site. An example of this type of docking software strategy is Wang et al[33]. [33] Wang, J,, Kolhnan, P. A. and Kuntz, I. D., "Flexible ligand docking: A multr step strategy approach", Proteins, Vol. 36,1-19 (1999). Force fields may be used to assign various atomic, bond, and/or other chemical or physical descriptors associated with components of molecules. In the context of the current invention this may include, but is not limited to, such items as vdW radii, charges (formal or partial), solvation dependent parameters, and equilibrium bond constants. An example cited in the teclmical description is the Tripos force field described in Clark et al. [34]. [34] Clark, M., Cramer, R.D., Opdenbosch, N. V., "Validation of the General Purpose Tripos 5.2 Force Field", J. Comp. Chem., Vol. 10,982-1012 (1989). BRIEF SUMMARY OF THE INVENTION Aspects of the present invention relate to a method and apparatus for the analysis of molecular combinations featuring two or more molecular subsets, wherein either one or both molecular subsets are from a plurality of molecular subsets selected from a molecule library, based on computations of shape complementarity utilizing a basis expansion representing ilecular shapes of the first and second molecular subsets in a coordinate system. Sets of asformed expansion coefficients are calculated for a sequence of different configurations, ,, relative positions and orientations, of the first molecular subset and the second molecular Dset using coordinate transformations. In order to first obtain sets of translated expansion efficients, the calculation includes application of a translation operator to a reference set of pansion coefficients before performing one or more rotation operations. The precomputed ts of translated expansion coefficients may be stored on a computer recordable medium, en later retrieved, and further subjected to one or more rotation operators in order to obtain te of transformed expansion coefficients. Then a shape complementarity score, presenting a correlation of the firet and second molecular subsets, is computed over the jquence of different sampled configurations for the molecular combination, where each ) impled configuration differs in both the relative positions and orientations of the first and 3cond molecular subsets. The aspect of the invention involving application of the translation perator prior to one or more rotation operator(s) has significant and beneficial implications or hardware-based implementations of the method. Various embodiments of the invention elating to efficient implementation in the context of a hardware apparatus are also discussed. PRIPF DESCRIPTION OF THE DRAWINGS A more complex appreciation of the invention and many of the advantages thereof mW be reachly obtained as the same becomes better understood by references to the detailed iescription when considered in cormection with the accompanying drawings, wherein: Fig. 1 is a block diagram view of an eatnbodiment of a system that utilizes the present invention in accordance with analysis of a molecular combinations based on computations of shape complementarity over a set of sampled configurations; Figs. 2a, 2b, and 2c respectively show a 'ball and stick' representation of a input pose for a methotrexate molecule, a digital representation in the form of a pdb formatted file, and another digital representation in the form of a mol2 formatted file, both files containing structural and chemical information for the molecule depicted in Fig. 2a; Fig. 3 shows a flow diagram of an exemplary method of assessing shape complementarity of two molecular subsets, performed in accordance with embodiments of the present invention; Fig. 4 shows illustrations of two molecular subsets assessed in accordance with embodiments of the present invention; Figs. 5 A and 5B. show two molecular subsets having internal and external volume ctions generated in accordance with embodiments of the present invention; Fig. 6 shows representations of two molecular subsets in discrete space for generating smal and external volume functions, in accordance with embodiments of the present 'ention; Fig. 7 illustrates how a 2-D continuous shape is discretized in accordance with ibodiments of the present invention; Fig. 8 shows coordinate-based representations of two molecular subsets in a joint ordinate system, in accordance with embodiments of the present invention; Fig. 9 shows the representation of various coordinate systems used in the present vention; Fig. 10 shows a representation of Euler angles as used in various embodiments of the resent invention. Fig. 11 shows a spherical sampling scheme used in embodiments of the present ivention. Fig. 12 is an illustration of spherical harmonics functions. Fig. 13 shows two molecular subsets in various configurations, i.e., having various elative translations and orientations to one another, for computing shape complementarity icores, in accordance with embodiments of the present invention; Fig. 14 shows a flow diagram of a method for analyzing a molecular combination x)mprising two molecular subsets based on computations of shape complementarity, by first performing a rotation operation on appropriate sets of expansion coefficients before applying a translation operator in order to generate sets of transformed expansion coefficients, as per the work of Ritchie et al; Fig. 15 shows a flow diagram of a novel, and substantially more efficient method, analyzing a molecular combination comprising two molecular subsets based on computations of shape complementarity, by first performing a translation operation on appropriate sets of expansion coeffficients before applying a rotation operator in order to generate sets of transformed expansion coefficients, in accordance with embodiments of the present invention. DETAILED DESCRIPTION OF THE INVENTION The present invention has many applications, as will be apparent after reachng this disclosure. In describing an embodiment of a computational system according to the present ivention, only a few of the possible variations are described. Other applications and ariations will be apparent to one of ordinary skill in the art, so the invention should not be jnstrued as narrowly as the examples, but rather in accordance with the appended claims. Embodiments of the invention will now be described, by way of example, not mitation. It is to be understood that the invention is of broad utility and may be used in lany different contexts. A molecular subset is a whole or parts of the components of a molecule, where the omponents can be single atoms or bonds, groups of atoms and/or bonds, amino add esidues, nucleotides, etc. A molecular subset might include a molecule, a part of a molecule, I chemical compound composed of one or more molecules (or other bio-reactive agents), a )rotein, one or more subsets or domains of a protein, a nucleic acid, one or more peptides, or me or more oligonucleotides. In another embodiment of the present invention, a molecular mbset may also include one or more ions, individual atoms, or whole or parts of other simple nolecules such as salts, gas molecules, water molecules, radicals, or even organic compounds ike alcohols, esters, ketones, simple sugars, etc. In yet another embodiment, the molecular subset may also include organic molecules, residues, nucleotides, carbohydrates, inorganic molecules, and other chemically active items including synthetic, medicinal, drug-like, or natural compounds. In yet another embodiment, the molecular subset may already be bound or attached to the target through one or more covalent bonds. In another embodiment the molecular subset may in fact include one or more structural components of the target, such as secondary structure elements that make-up a tertiary structure of a protein or subunits of a protein quaternary structure. In another embodiment the molecular subset may include one or more portions of a target molecule, such as protein domains that include the whole or part of an active site, one or more spatially cormected subsets of the protein structure that are selected based on proximity to one or more protein residues, or even discormected protein subsets that feature catalytic or other surface residues that are of interest for various molecular interactions. In another embodiment, the molecular subset may include the whole of or part of an existing molecular complex, meaning a molecular combination between two or more other molecular subset, as, for example, an activated protein or an allosterically bound protein. A molecular combination (or combination) is a collection of two or more molecular subsets that may potentially bind, form a molecular complex, or otherwise interact with one nother. A combination specifies at the very least the identities of the two or more iteracting molecular subsets, A molecular pose is the geometric state of a molecular subset described by its position nd orientation within the context of a prescribed coordinate system. A molecular onfiguration (or configuration) of a molecular combination represents the joint poses of all onstituent molecular subsets of a molecular combination. Different configurations are lenoted by different relative positions and orientations of the molecular subsets with respect 0 one another. Linear coordinate transformations that do not change the relative position or mentation of constituent molecular subsets will not result in different configurations. For the purposes of the invention different configurations of a molecular combination ire obtained by the application of rigid body transformations, including rdative translation md rotation, to one or more molecular subsets. For the purposes of the invention, such rigid body transformations are expected to preserve the conformational structure, as well as the stereochemistry and/or tautomwism (if applicable), of each molecular subset, fa regards to the invention it is contemplated that when analyzing distinct conformations or stereoisomers of a molecular subset, each distinct conformation or stereoisomer will appear in a distinct molecular combination, each with its own attendant analysis. In this way, molecular combinations featuring flexible molecular subsets may be better analyzed using the invention based on consideration of multiple combinations comprising distinct conformations and/or stereoisomers. In many of the forthcoming examples and explanations, the molecular combination will represent the typical scenario of two molecular subsets where a ligand biomolecule (first molecular subset) interacts with a target biomolecule (usually a biopolymer; second molecular subset). Thus in regards to the present invention, an of a molecular combination may seek to determine whether, in what fashion (i.e., binding mode), and/or to what degree, a ligand will interact with a target molecule based on computations of shape complementarity of one or more configurations. Here the term "shape complementarity" measures the geometric fit or correlation between the molecular shapes of two molecules. The concept can be generalized to any two objects. For example, two pieces of a jigsaw puzzle that fit each other exhibit strong shape complementarity. It should be understood that, unless otherwise indicated, such examples and explanations could more generally apply to molecular combinations wherein more than two molecular subsets bind or interact with one another, representing the whole of, or portion(s) of, one or more target molecules and/or one or more ligands. As an example, in one embodiment of the present invention the molecular combination may represent a target interacting with a ligand (i.e., target-ligand pair) where )ne molecular subset is from the protein and the other the ligand. In a further embodiment, he molecular combination may represent a target-ligand pair where one molecular subset is he entire ligand biomolecule but the other molecular subset is a portion of a target biopolymer containing one or more relevant active sites. In yet another embodiment, the molecular combination may feature more than two nolecular subsets, one representing a target (whole or part) and the other two correspond to two distinct ligands interacting with the same target at the same time, such as in the case of competitive thermodynamic equilibrium between a possible inhibitor and a natural binder of a protein. In yet another embodiment the previous example may be turned around such that the molecular combination features two target molecules in competition with one ligand biomolecule. As another example, in one embodiment the molecular combination may represent a proteiirrprotein interaction in which there are two molecular subsets, each representing the whole or a relevant portion of one protein. In a further embodiment, the molecular combinations may also represent a proteirrprotein interaction, but now with potentially more than two molecular subsets, each representing an appropriate protein domain. As a fiffther example, the molecular combination may feature two molecular subsets representing a target-ligand pair but also additional molecular subsets representing other atoms or molecules (heteroatoms or heteromolecules) relevant to the interaction, such as, but not limited to, one or more catalytic or structural metal ions, one or more ordered, bound, or structural water molecules, one or more salt molecules, or even other molecules such as various lipids, carbohydrates, acids, bases, mRNA, ATP/ADP, etc. In yet another embodiment, the molecular combination may feature two molecular subsets representing a target-ligand pair but also one or more added molecular subsets representing a whole or portion of a cell membrane, such as a section of a lipid brlayer, nuclear membrane, etc., or a whole or portion of an organelle such as a mitochondrion, a ribosome, endoplasmic reticulum, etc. In another embodiment, the molecular combination may feature two or more molecular subsets, with one or more molecular subsets representing various portions of a molecular complex and another subset representing the ligand interacting with the complex at an unoccupied active site, such as for proteins complexed with an allosteric activator or for proteins containing multiple, distinct active sites. ' In another embodiment, the molecular combination may feature two or more molecular subsets representing protein chains or subunits interacting noncovalently as per a uatemary protein structure. In another embodiment, the molecular combination may feature vo or more molecular subsets representing protein secondary structure elements interacting s per a tertiary structure of a polypeptide chain, induced for example by protein folding or mutagenesis. In many of the forthcoming examples and explanations, the molecular combination /ill represent the typical scenario of a target-ligand pair interacting with one another. As Iready mentioned in regards to the present invention, an analysis of a molecular combination nay seek to determine whether, in what fashion, and/or to what degree or with what ikelihood, a ligand will interact with a target molecule based on computations of shape ;omplementarity. In another embodiment, the analysis may involve a plurality of molecular combinations, each conresponding to a different ligand, selected, for example, from a nolecule library (virtual or otherwise), in combination with the same target molecule, in Drder to find one or more ligands that demonstrate high shape complementarity with the target, and are therefore likely to bind or otherwise react with the target. In such cases, it tnay be necessary to assign a score or ranking to each analyzed molecular combination based on the estimated maximal shape complementarity across a set of different configurations for eadh combination, in order to achieve relative comparison of relevant predicted bioactivity. In such a scenario where each target-ligand pair is an individual combination, and if there areiV^Hgands to be tested against one target, then there will be A^ distinct molecular combinations involved in the analysis. For sufficiently large molecule libraries, it may be necessary to analyze millions or more potential molecular combinations for a single target protein. In yet another embodiment, the analysis may be reversed and the plurality of molecular combinations represents a plurality of target molecules, each in combination with the same ligand biomolecule in the same environment. In other embodiments, the molecular combinations may represent multiple ligands and/or targets reacting simultaneously, i.e., more than just a target-ligand pair, and may also include various heteroatoms or molecules as previously discussed. Fig. 1 illustrates a modeling system 100 for the analysis of molecular combinations based on computations of shape complementarity across a set of configurations for the molecular combination. As shown a configuration analyzer 102 receives one or more input (or reference) configuration records 106, including relevant structural, chemical, and physical data associated with input structures for both molecular subsets from an input molecular combination database 104. The configuration analyzer 102 comprises a configuration data ansformation engine 108 and a shape complementarity engine 109. Results fi'om the Dnfiguration analyzer 102 are output as configuration results records 111 to a configuration jsults database 110. Modeling system 100 may be used to efficiently analyze molecular combinations via omputations of shape complanementarity. In some embodiments, this may include, but is not imited to, prediction of likelihood of formation of a potential molecular complex, or a proxy liereof, the estimation of the binding affnity between molecular subsets in the molecular ombination, the prediction of the binding mode (or even additional alternative modes) for he molecular combination, or the rank prioritization of a collection of molecular subsets e.g., ligands) based on maximal shape complementarity with a target molecular subset across iampled configurations of the combination, and would therefore also include usage issociated with computational target-ligand docking. Furthermore, the method provides for performing a dense search in the uonfigurational space of two or more molecular subsets having rigid bothes, that is, assessing relative orientations and translations of the constituent molecular subsets. The method can also be used in conjunction with a process for generating likely yet distinct conformations of one or both molecular subsets, in order to better analyze those molecular combinations where one or both of the molecular subsets are flexible. In a typical operation, many molecular combinations, each featuring many different configurations, may be analyzed. Since the total possible number of configurations may be enormous, the modeling system 100 may sample a subset of configurations during the analysis procedure according to an appropriate sampling scheme as will be discussed later. However, the sampled subset may still be very large (e.'g., millions or even possibly billions of configurations per combination). A shape complementarity score is generated for each sampled configuration and the results for one or more configurations recorded in a storage medium. The molecular combination may then be assessed by examination of the set of configuration results including the corresponding computed shape complementarity scores. Once the cycle of computation is complete for one molecular combination, modeling of the next molecular combination may ensue. Alternatively, in some embodiments of the modeling system 100, multiple molecular combinations may be modeled in parallel. Likewise, in some embodiments, during modeling of a molecular combination, more than one configuration may be processed in parallel as opposed to simply in sequence. In one embodiment, modeling system 100 may be implemented on a dedicated licroprocessor, ASIC, or FPGA. In another embodiment, modeling system 100 may be nplemented on an electronic or system board featuring multiple microprocessors, ASICs, or PGAs. In yet another embodiment, modeling system 100 may be implemented on or across wltiple boards housed in one or more electronic devices. In yet another embodiment, lodeling system 100 may be implemented across multiple devices containing one or more licroprocessors, ASICs, or FPGAs on one or more electronic boards and the devices ormected across a network. In some embodiments, modeling system 100 may also include one or more storage nedia devices for the storage of various, required data elements used in or produced by the inalysis. Alternatively, in some other embodiments, some or all of the storage media devices I nay be externally located but networked or otherwise cormected to the modeling system 100. Examples of external storage media devices may include one or more database servers or file systems. In some embodiments involving implementations featuring one or more boards, the modeling system 100 may also include one or more software processing components in order to assist the computational process. Alternatively, in some other embodiments, some or all of the software processing components may be externally located but networked or otherwise cormected to the modeling system 100. hi some embodiments, results records from database 110 may be further subjected to a configuration selector 112 during which one or more configurations may be selected based on various results criteria and then resubmitted to the configuration analyzer 102 (possibly under different operational conditions) for further scrutiny (i.e., a feedback cycle). In such embodiments, the molecular configurations are transmitted as inputs to the configuration analyzer 102 in the form of selected configuration records 114. In another embodiment, the configuration selector 112 may examine the results records from database 110 and construct other configurations to be subsequently modeled by configuration analyzer 102. For example, if the configuration analyzer modeled ten target-ligand configurations for a given target-ligand pair and two of the configurations had substantially higiher estimated shape complementarity than the other eight, then the configuration selector 112 may generate further additional configurations that are highly similar to the top two high-scoring configurations and then schedule the new configurations for processing by configuration analyzer 102. In some embodiments, once analysis of a molecular combination is completed (i.e., all desired configurations assessed) a combination post-processor 116 may used to select one or lore configuration results records from database 110 in'order to generate one or more either ualitative or quantitative measures for the combination, such as a combination score, a Dmbination summary, a combination grade, etc., and the resultant combination measures are len stored in a combination results database 118. In one embodiment, the combination leasure may reflect the configuration record stored in database 110 with the best-observed bape complementarity. In another embodiment, multiple configurations with high shape omplementarity are submitted to the combination post-processor 116 and a set of ombination measures written to the combination results database 118. In another mbodiment, the selection of multiple configurations for use by the combination post->rocessor 116 may involved one or more thresholds or other decisiorrbased criteria. In a further embodiment, the combination measures output to the combination results latabase 118 are based on various statistical analysis of a sampling of possibly a large ramber of configuration results records stored in database 110. In other embodiment the selection sampling itself may be based on statistical methods (e.g., principal componoit malysis, multidimensional clustraing, multivariate regression, etc.) or on patterrrmatching methods (e.g., neural networks, support vector machines, etc.) In another embodimeot, the combination post-processor 116 may be applied dynamically (i.e., orrthe-fly) to the configuration results database 110 in parallel with the analysis of the molecular combination as configuration results records become available. In yet another embodiment, the combination post-processor 116 may be used to rank different configurations in order to store a sorted list of either.all or a subset of the configurations stored in database 110 that are associated with the combination in question. In yet other embodiments, once the final combination results records, reflecting the complete analysis of the molecular combination by the configuration analyzer 102, have been stored in database 118, some or all of the configuration records in database 110 may be removed or deleted in order to conserve storage in the context of a library screen involving possibly many different molecular combinations. Alternatively, some form of garbage collection may be used in other embodiments to dynamically remove poor configuration results records from database 110. In one embodiment, the molecular combination record database 104 may comprise one or more molecule records databases (e.g., flat file, relational, object oriented, etc.) or file systems and the configuration analyzer 102 receives an input molecule record corresponding to an input structure for each molecular subset of the combination. In another embodiment, when modeling target proteirrligand molecular combinations, the molecular combination scord database 104 is replaced by an input target record database and an input ligand (or rug candidate) record database. In a further embodiment, the input target molecular records lay be based on either experimentally derived (e.g., X-ray crystallography, NMR, etc.), nergy minimized, or model-built 3-D protein structures. In another embodiment, the input gaud molecular records may reflect energy minimized or randomized 3-D structures or other -D structures converted from a 2-D chemical representation, or even a sampling of low nergy conformers of the ligand in isolation. In yet another embodiment, the input ligand lolecular records may correspond to naturally existing compounds or even to virtually ;enerated compounds, which may or may not be synthesizable. In order to better illustrate an example of an input structure and the associated input nolecule record(s) that may form an input configuration record submitted to configuration inalyzer 102 we refer the reader to Figs. 2a, 2b, and 2c. Fig. 2a shows a "ball-and-stick" rendering of a pose 205 of a methotrexate molecule iOO with chemical formula CaoHaaNgOs. The depicted molecular subset comprises a joUection of atoms 220 and bonds 230. The small, black atoms, as indicated by item 213, represent carbon atoms. The tiny, white atoms, as indicated by item 216, represent hydrogen atoms, whereas the slightly larger dark atoms (item 210) are oxygen atoms and the larger white atoms (item 229) are nitrogen atoms. Continuing in Fig. 2a, item 223 denotes a circle containing a benzene ring (C6H4), and item 225 a circle containing a carboxyl group (COO"), and itism 227 another circle containing a methyl group (CH3). Item 233 denotes a covalent bond cormecting the benzene ring 223 to the ester group that includes the methyl group 227. Item 235 denotes a covalent bond cormecting the carbon atom 213 to the carboxyl group 225. Lastly item 237 denotes a covalent bond cormecting the methyl group 227 to a nitrogen atom 229. Fig. 2b shows a pdb file representation 240 of a chemical structure for liie methotrexate ligand pose described in Fig. 2a, including a general header 250, a section 260 composed of atom type and coordinate information, and a section 265 regarding bond coimectivity information. The header section 250 may contain any armotation or other information desired regarding the identity, source, or characteristics of the molecular subset and its conformation and/or stereochemistry^ Section 0260 shows a list of all 33 nonhydrogen atoms of methotrexate and for each atom it includes a chemical type (e.g., atomic element) and three spatial coordinates. For instance, the line for atom 6 shows that it is a nitrogen atom with name NA4 in a compound (or residue if a protein) named MTX in chain A with compound (or residue) id of 1 and with (x, y, z) coordinates (20.821,57.440, .075) in a specified Cartesian coordinate system. Note that the compound or residue name sld may be more relevant for amino or nucleic acid residues in biopolymers. Section 265 of the PDB file 240, sometimes called the cormect record of a PDB file, sscribes a list of the bonds associated with each atom. For instance, the first line of this iction shows that atom 1 is bonded to atoms (2), and (12), whereas the second line shows vat atom 2 is bonded to atoms (1), (3), and (4). Notice also how in this example hydrogens re missing and as such the bond cormections for each atom may not be complete. Of course, ompleted variants of the PDB file representation are possible if the positions of hydrogen toms are already specified, but in many cases where the chemical structure originates from sxperimental observations the positions of hydrogens may be very uncertain or missing iltogether. Fig. 2c shows a Tripos mol2 file containing various additional chemical descriptors ibove and beyond the information shown in the PDB file in Fig. 2b. Column 270 lists an index for each atom; column 273 lists an atom name (may be nonunique) for each atom; columns 275,277, and 279 respectively list x, y, z coordinates for each atom in an internal coordinate system; column 280 lists a SYBYL atom type according to the Tripos force field [34] for each atom that codifies information for hybridization states, chemical type, bond cormectivity, hydrogen bond capacity, aromaticity, and in some cases chemical group; and columns 282 and 285 list a residue id and a residue name for each atom (relevant for proteins, nucleic acids, etc.). Section 290 lists all bonds in the molecular subset. Column 291 lists a bond index for each bond; columns 292 and 293 the atom indices of the two atoms cormected by the bond; and column 295 the bond type, which may be single, double, triple, delocalized, amide, aromatic, or other specialized covalent bonds. In other embodiments such information may also represent noncovalent bonds such as salt bridges or hydrogen bonds. In this example, notice how the hydrogen atoms have now been included. In one embodiment the configuration data transformation engine 108 may directly transform one or more input molecular configurations into one or more other new configurations by application of various rigid body transformations. In other embodiments, the configuration data transformation engine 108 may instead apply rigid body transformations to sets of basis expansion coefficients representing molecular shapes for reference poses for each molecular subset as will be discussed in more detail later in the technical description. In some embodiments, the set of configurations visited during the course of an analysis of a molecular combination may be determined according to a schedule sampling scheme specified in accordance with a search of the permitted configuration )ace for the molecular combination. In some embodiments, whether generated by direct transformation of structural jordinates or by transformation of sets of basis expansion coefficients, the configuration ata transformation engine 108 may produce new configurations (or new sets of basis Kpansion coefficients corresponding to new configurations) sequentially and feed them to le shape complementarity engine 109 in a sequential marmer, or may instead produce them 1 parallel and submit them in parallel to the shape complementarity engine 109. The shape complementarity engine 109 is responsible for gmerating a shape iomplementarity score or equivalent measure for each sampled configuration of the nolecular combinations and makes use of the present invention to efficiently compute the ;hape complementarity for each configuration based on use of basis expansions and rigid jody transformations of molecular shapes. The shape complementarity engine 109 may also include one or more storage components for data specific to the computations of shape complementarity. In some embodiments,, the configuration results records 111 may include a quantitative measure related to the shape complementerty evaluated for each configuration. In one embodiment, this may be a score. In another embodiment, this may be a probability. In other embodiments, the configuration results records 111 may include a qualitative measure related to the shape complementarity evaluated for the configuration. In one embodiment, this may be a grade. In another embodiment this may be a categorization (i.e., poor, weak, strong, etc.). In yet another embodiment this may be a simple pass-fail measure. In many embodiments, the configuration results records 111 may also include information used to specify the identity and/or nature of configuration corresponding to a given shape complementarity score. In addition to the identity of the interacting molecular subsets, there may be a need to armotate or otherwise represent the geometrical state of the configuration. Typically this may be achieved by storing the parameters of the rigid body transformation used to generate the configuration from an input or reference configuration. In some embodiments, the configuration selector 112 may utilize various selection criteria in order to resubmit certain configurations back to modeling system 102 for more computations. In one embodiment, the selection criteria may be predicated on passing of a threshold or other decision mechanism based on one or more qualitative affinity measures. In other embodiment, the selection criteria may be based on a threshold or other decision jchanism based on one or more quantitative shape complementarity scores. In yet another embodiment, the selection criteria used by the configuration selector .2 may be based on various statistical analysis of a number of different configuration results cords stored in database 110, including, but not limited to, principal component analysis, lultidimensional clustering, Bayesian filters, multivariate regression analysis, etc. In yet lother embodiment, the selection criteria may be based on various pattern matching analysis fa number of different configuration results records stored in database 110, including, but ot limited to, use of neural networics, support vector machines, hidden Markov models, etc. In some embodiments, the configuration data transformation engine 108 may receive sertain resubmitted configurations from the configuration selector 112 and utilize them as nputs to start a new cycle of shape complementarity computations. For example, il a particular configuration was selected from database 110 based on high shape complementarity by the configurations selector 112, the configuration data transformation engine 108 may generate multiple configurations (or multiple sets of basis expansion coefficients corresponding to new configurations) that are similar (i.e., slightly dififwwit positions and orientations for each molecular subset) in order to better investigate that portion of the possible configuration space of the molecular combination. In other embodiments, the new cycle of shape complementarity computations instigated by fee resubmission of the selected configurations records 114 may involve the operation of the configuration analyzer ' 102 under a different set of conditions or using a different set of control parameters. In further «nbodiments, the selected configurations records ljl4 may kick off a new cycle of shape complementarity computations using a different variant of the configuration analyzer 102, including the use of a modified formulation for subsequent shape complementarity scores (if appropriate). Fig. 3 shows a flow diagram of an exemplary method 300 of analyzing a molecular combination based on computations of shape complementarity across a set of configurations, performed in accordance with embodiments of fee present invention. The mefeod 300 of Fig. 3 is described wife reference to Figs. 4-15. As explained below, fee method 300 generally involves computing a basis expansion representing molecular shapes of fee constituent molecular subsets, computing transformed expansion coefficients for different configurations (i.e., relative positions and orientations) of fee molecular subsets, and computing a correlation function representing a shape complementarity of the two molecular ibsets using the transformed expansion coefficients. Embodiments of this method corporate various combinations of hardware, software, and firmware to perform the steps sscribed below. In Fig, 3, in step 302, a first molecular subset 410 and a second molecular subset 420 re provided, as shown in Fig. 4. Each molecular subset has a molecular shape, illustrated in ig. 4, As used herein, "molecular shape" generally refers to a volumetric function spresenting the structure of a molecular subset comprising a plurality of atoms and bonds. Tiose skilled in the art will appreciate that the molecular subsets may have various shapes itW tiian those shown in Fig. 4. The first molecular subset 410 has a plurality of atomis and )onds. Some of these atoms in the first molecular subset are surface atoms 425. The surface itoms are proximal to and define a molecular surface 430 for the first molecular subset 410 jased on the location ofthose surfaced atoms 425, Similarly, the second molecular subset 420, as shown in Fig, 4, also has a plurality of atoms and bonds. Some of the atoms in the second molecular subset are surface atoms 435, the locations of which define a molecular surface 440 for the second molecular subset 420. The molecular surface 430 of the first molecular subset can be a solvwit accessible molecular surface, which is generally the surface traced by the center of a small sphere rolling over the molecular surface 430. As used herein, "solvent" gaierally refers to the plurality of atoms, ions, and/or simple molecules (e.g., water, salt, sugars) that comprise an ambient medium, polarizable or otherwise. Computational methods for generation of solvent accessible surfaces includes the method presented in Cormolly, M. L., "Analytical molecular surface calculation.", (1983), J. Applied Crystallography. 16,548-558; all of which is hereby incorporated in its entirety. In Fig. 4, often the first molecular subset 410 is a ligand, and the second molecular subset 420 is a protein. However, as already discussed in regards to the definition of a molecular subset, the molecular subsets 410 and 420 can have various compositions. In Fig. 3, in step 304, a first internal volume function is generated. This first internal volume function is a representation of the subset of a volume enclosed by the first molecular surface 430 of molecular subset 410. As used herein, the "internal volume function" of a molecular subset is graierally any subset of the 3-D volume enclosed by the molecular surface. In one embodiment, the first intemal volume function is defined as a union of a set of kernel functions, where generally a kernel function is a 3-D volumetric function with finite support in a localized neighborhood about an atom and/or bond, and where each kernel function is associated with atoms and/or bonds in the first molecular subset 410. In one embodiment, the kernel function used in defining an internal volume function a 3-D Gaussian function localized around the center of an atom. In other embodiments, ich kernel function can be dependent on the chemical identity of associated atoms and/or Dnds. Alternatively, the kernel function can be dependent on the location of the associated torn or bond within a themical group. In one embodiment, the kernel function associated rith an atom is a nonzwo constant for positions within a Van der Waals (VdW) sphere, i.e., a phere with radius equal to the VdW radius of an atom of the given type and centered on the torn, and has a value of zero at other positions. This nonzoro constant can have a value of inity, for example. In another embodiment, the first internal volume function is further specified by ntersection with a proximal surface volume defined by the movement of a probe sphere 520 white balls), as shown in Fig. 5A, which movers at positions proximal to and internal to a irst molecular surface 510 for the first molecular subset 410. The proximal surface volume IS shown in Fig. 5B as the shaded region 540, which depicts the same molecular subset 410, and fee first internal volume function resulting from the intersection of the molecular subset and region 540 is displayed as the volume occupied by the dark balls 550 also in Fig. 5B. In yet another embodiment, the first intern^ volume function is wholly comprised of the proximal surface volume depicted in Fig. 5B as the shaded region 540. Similarly, a second internal volume function may be generated for the second molecular subset 420 and can be defined in a similar marmer as the first internal volume function. Also in step 304 of Fig. 3, external volume functionss are defined for both the first molecular subset and the second molecular subset. As used h«ran, the "external volume function" of a molecular subset is generally any subset of the 3-D volume external to the molecular surface of that molecular subset. A first extern^ volume function is generated as a representation of a subset of a volume external to the first molecular surface. As shown in Fig. 5A, this extemal volume function can be defined by the movement of a probe sphere 530 (white balls) that moves at positions proximal to and extarfial to the first molecular surface 510. Similarly, a second extemal volume function is defined for the second molecular subset In one embodiment, when the volume functions are generated, the probe spheres 520 (internal) or 530 (extemal) travels along the entire molecular surface 510, while in another embodiment, the probe sphere moves along a portion of the molecular surface 510. In one embodiment, the probe sphere has a constant radius, while in another embodiment the radius f the probe sphere varies as a function of the location of the probe sphere on the molecular jrface. In Fig. 5B, the external volume function and the internal volume function for first lolecular subset 410 are shown. In Fig. SB, the region 550, represented as a union of dark alls, identifies the volumetric domain of the internal volume function, x, while the diagonal ashed region 560 represents the volumetric domain of the external volume function, a, of le first molecular subset 410. The volumetric domains for T and a may be similarly ientified for the second molecular subset 420. The volumetric functions, T and a, are used D calculate shape correlations, as described below. In one embodiment, described with reference to Fig. 6, the first and second molecular ubsets 410 and 420 are represented in discrete space 600 in order to generate the respective ntemal and extemal volume functions. As used herein, "discretization" generally refers to sonverting a continuous representation to a discrete one, e.g., converting the function from its continuous representation into a series of numbers that best approximates the continuous iunction as projected onto a set of grid cells. For example, the Fig. 7 illustrates how a general 2-D continuous shape is discretized in a 2-D rectilinear grid. The black dots represent centers of the occupied grid cells; the ivhite dots represeint cwiters of unoccupied grid cells. In Fig. 6, the extemal volume function, a, for a particular grid cell is assigned a nonzero numerical value when a grid cell is inside the volumetric domain 660 of Fig. 6, i.e., the grid cell is occupied, and zero otherwise. Similarly, the internal volume function, T, for a particular grid cell is assigned a nonzero numerical value when the grid cell is inside the volumetric domain 650 of Fig. 6, and zero otherwise. In one embodiment only a significant fraction of the grid cell must lie within the appropriate domain in order for the grid cell to be considered occupied. While Fig. 6 shows a two-dimensional cross-sectional view of the volume functions for the molecular subsets and a 2-D Cartesian grid, those skilled in the art should understand that the principles described above are equally applicable to three-dunensional and higher multidimensional spaces, as well as to other coordinate based representations, where the phrase "coordinate based representation" generally refers to representing a function in terms of coordinates of a coordinate system. In one embodiment, a Cartesian coordinate based representation is used where each grid cell in three-dimensions is a cuboid. The cuboid grid cells in Fig. 6 with nonzero values •r cr are illustrated as having horizontal stripes, and those with nonzero values for t are lustrated as having vertical stripes. In this way, values are assigned to a and T, for each grid 5II so that the internal and external volume functions are represented as a set of numbers for le entirety of grid cells. In one embodiment, the positive numerical value assigned when a rid cellis occupied is unity and thus the representative set of numbers is a 3-D bitmap. In Fig. 3, in step 306, the method 300 includes defining individual coordinate-based epresentations for the molecular subsets 410 and 420, wherdiy each molecular subset is epresented in a coordinate system. A three-dimensional coordinate system is a systematic vay of describing points in three-dimensional space using sets of three numbers (or points in I plane using pairs of numbCTS for a two-dimensional space). As shown in Fig. 8, an ndi vidual coordinate based representation of the first molecular subset 410 is defined using a iirst coordinate system 810. An individual coordinate based representation of the second molecular subset 420 is defined using a second coordinate system 820. In one embodiment, the individual coordinate based representations are defined using a spherical polar coordinate system. The spherical polar coordinate system is a three-dimensional coordinate system where the coordinates are as follows: a distance from the origin r, and two angjes 9 and cp found by drawing a line from the given point to the origin and measuring the angles formed with a given plane and a givra line in that plane. Angle 6 is taken as the polar (co-latitudinal) coordinate with 9 €[0, n] and angle cp is the azimuflial (longitudinal) coordinate with

= 30. In anotha: embodiment, the order of the expansion is adaptively determined based on a preliminary quantitative analysis of represKitation errors for trial values of the expansion order, and may therefore be of differwit magnitude for different pairs, of molecular subsets 410 and 420 based on the characteristics of their respective internal and external volume functions. In one embodiment, the basis expansion is an orthogonal basis expansion comprising a plurality of mutually orthogonal basis functions. If the basis functions satisfy the following mathematical condition, they are called mutually orthogonal: here Cjj is a constant (not necsessarily unity when i = j), 8ij is the usual Kronecker delta, and ,e integral is over the entire M-dimensional space. For an orthogonal basis expansion, an expansion coefficient, % corresponding to a articular basis function, B,, can be written as follows: vhere Cii is a constant. However, once again for the practical pxirposes of computation, the expansion x>e£fici«its are discaretized by converting the integral in Eqn. 7 to a finite summation. In the 3ase of a set of expansion coefiScients for an orthogonal basis expansion, the discretized expansion coefficient, aj, for an orthonormal basis function, Bi, takes the following form: where the siunmation is over the discrete points c, i.e., x^ is a sample point in the M- dimensional space represented here by x. In another embodiment, the basis expansion is an orthonormal basis expaiision comprising a plurality of mutually orthonormal basis functions. If the basis functions are mutually orthogonal and in Eqn. 8, Qi is unity for all relevant basis functions, then the basis functions are said to be mutually orthonormal. This similarly simplifies the expressions for ai in eqns. 7 and 8. A general 3-D function in spherical polar coordinates can be represented in terms of a radial/spherical harmonics basis expansion comprising a plurality of basis functions, each basis function defined as the product of one of a set of orthonormal radial basis functions, R^(r), and one of a set of real-valued spherical harmonics basis functions, y!"(6,(|>), as follows: here {anim} is the set of radial / spherical harmonics expansion coefficients, (r, 9, ^) are the iherical coordinates of a point in 3D space, n = [1, N), integer, I = [0, rr1], integer, m = [-1, , integer. The usage of such an expansion is common practice in the quantum mechanical iscription of numerous atomic and molecular orbitals. Hence the indices n, 1, and |m| S: 0 •e often respectively referred to as the principal quantum numbor, angular quantum (or rbital) number and azimuthal (or magnetic moment) quantum number. In Eqn. 9, each radial basis function, R^(r), is a 1 -D orthonormal basis function epending solely on the radius, r. , The form for the radial basis functions is often chosen based on the problem at hand, .g., the scaled hydrogen atom radial wave function in quantum mechanics is based for ixample on associated Laguerre polynomials (Arfken e? a/.) as follows: ivhere the square root term in the normalization factor, p is the scaled distance, p = i^/k, k is the scaling parameter, F is the gamma function, and LQ are the associated Lagu^re polynomials; where a general Laguerre polynomial is a solution to the Laguerre differantial equation given by: xy" + (1 - x)y' + Ay = 0 [Eqn. ii] and the associated Laguerre polynomials themselves are given explicitly by Rodrigue's formula as follows: Various radial basis functions can be used in accordance with embodiments of the present invention. In one embodiment, the radial basis functions include fee scaled Laguerre polynomial-based functions of Eqn. 10. In another embodiment, the radial basis functions include unsealed forms of Eqn. 10 in terms of r (not p) and without the normalization »nstants. In yet another embodiment the radial basis functions include a Bessel function of e first kind (Jn(r)). In yet another embodiment the radial basis functions include a Hermite slynomial function (Hn(r)), In other embodiments, the radial basis functions include any lutually orthonormal set of basis functions that depend on the radius in a spherical ^ordinate system centered on the respective molecular center of the molecular subset in uestion. In Eqn. 9, each real-valued spherical harmonic basis function, y!"(0,<^), is a 2-D rthonormal basis function depending on the angular variables (6, ), can be obtained from suitable linear combinations of JT and its complex conjugate Y*i" in ordtt" to represent the real and imaginary parts of ll*", as follows: ased on Eqn. 9, the expansion coefficients {anim} coefficients for an arbitrary 3-D volume motion f(r, 0, <|)) are defined as: /here the integral is over the extent of functions f(r, 6, (|») in spherical coordinates and dV is differential volume element in spherical coordinates. The discretized analog of the expansion coefficients in Eqn. 17 are given by: (vhexe the summation is over all grid cells, c, and (r^, 9^, ^^ are the spherical coordinates of the center of grid cell c and AVc is the volume of grid cell c. In one embodiment of the present invention, the gad cells in Eqn, 18 are cuboids from a Cartesian coordinate system and the spherical 3-tuples (r^, 0g, ^^ are converted 'orrtibe-fly' to Cartesian 3-tuples (Xc, yc, Zc) by means of a suitable coordinate transformation, for easy addressing of the function, f, over a lattice representation stored in a compute readable manory. In another embodiment the grid cells can be a varying volume, AVc. In yet another embodiment thp grid cells represent small volumes in a spherical coordinate system. In yet another embodiment the grid cells represent small volumes in 8 cylindrical coordinate system. Eqn. 18 has been used to represent volumetric functions desoibing the shape of a molecular subset by a corre^onding set of expansion coefficients as described in Ritchie et ah In one embodiment of the present invrention, in step 316 of Fig. 3, Eqn. 18 is used to represent the coordinate based representation of the intepnal volume function of the first molecular subset 410, ta, in temis of a set of expansion coefficients for an initial pose of molecular subset 410, herein designated as a reference set of expansion coefficients. Similarly, the coordinate based representation of the external volume function of the first molecular subset 410, CTO. is represented in terms of a set of expansion coefficients for the same initial pose of molecular subset 410, also designated as a refwence set of expansion coefficients but now for the external volume function, Og. In another embodimeait, fCr^, 6^, ^^ = 1 if the grid cell is occupied, i.e., "lies within" (see above) the nonzero domain of the volume function in question, and zero otherwise. Similarly, reference sets of expansion coefficients can be constructed for the internal d external volume functions of molecular subset 420 (respectively x\, and Ob) for an initial (se of molecular subset 420. Thus altogether there are four sets of reference expansion ^efficients computed for the two molecular subsets 410 and 420, each corresponding to a )lume function for one of the molecular subsets. The coefficient sets for molecular subset | 10 are designated as {a'nim} and {a^„im}, respectively, and the coefficient sets for molecular ibset 420 are designated as {h^tom} and {b^im} respectively. In another embodiment, the )ur reference sets of expansion coefficients are computed using Eqn. 8 where the set of basis motions {B;} correspond to a general basis .expansion, i.e., need not be the radial/sphaical armonics expansion of Eqn. 9, upon which Eqn. 18 is predicated. In one embodiment of the present invaition, the set of values comprising {f(jr^,9^, ig)} for all grid cells c (both occupied and unoccupied) in a coordinate based represrentation or the internal volume function of molecular subset 410 are converted to a stream or an array )f Cartesiarrbased values {f (x^, y^, z^)} via a suitable coordinate transformation and stored )n a computer readable and recordable medium for future retrieval. Later, when the stored /alues are to be used in the context of Eqn. 18 to compute expansion coefficients, the stored values are first retrieved and then converted back into {^r^, 0^, ^J) by an inverse coordinate transformation. In the case that fit^, 9^, i^^ is unity for occupied grid cells and zero otherwise, the stored set of values {f (x^, y^,, z^)} become a bit stream or a bitmap. Similarly the values corresponding to coordinate based representations for the external volume function of molecular subset 410 can be stored and retrieved in a similar marmer. The same also applies to the values corresponding to coordinate based representations for the internal and external volume functions of molecular subset 420. In Fig. 3, in step 318, the method continues wife providing a translation operator representing translation of the coordinate based representation 810 of the &st molecular subset with respect to the coordinate based representation 820 of the second molecular subset in the joint coordinate system. The term "translation opaator" refers to an operator that, when applied to a point, results in the point's translation along a vector as defined by the translation operator. The operator can be applied to any collections of points as a subset of 3-D space, e.g., a line, a curve, a surface, or a volume. In one embodiment, fee translation operator is a matrix function of fee displacement along fee intermolecular axis 870 between fee first and second molecular subsets. Then the anslation operator can be directly applied to the set of reference expansion coefficients for le internal and external volume functions of molecular subset 410, while leaving molecular jbset 420 untouched. In another embodiment, molecular subset 410 is held fixed and the anslation operator is directly applied to the set of reference expansion coefficiaits for the rtemal and external volume functions of molecular subset 420. In one embodiment, using the joint (R,Pi,Yi.o'2,p2,Y2) coordinate system of Fig. 8, in onjunction with the radial/spherical harmonics expansion of Eqn. 9, the translation operator epresenting the translation of the coordinate based representation of the internal volume unction for molecular subset 410 from R=0 (meaning molecular centers of molecular subsets fl 0 and 420 are same point in Fig. 8) to R > 0 is directly applied to the reference set of ixpansion coefficients {a^'^nim} for the internal and extfflnal volume functions for molecular ubset 410 according to the following rule: pvhere jSpf^ j are the new translated set of expansion coefScients, (n, 1, m) are quantum lumbers for the new translated expansion coefficient, (n', 1', m') are quantum numbers for the :>ne of the set of reference expansion coefficients, fht sununation is over all possible values of n' and 1', 5m„,- is the standard Kronecker delta, and K^-n'imi are matrix elements of a translation matrix function with values equal to resultant overly integrals betweai two different basis functions of the radial / spherical harmonics expansion, with quantum numbers n, 1, m and n', 1', m' respectively, separated by a distance R, and which are nonzero only when m = m'. The exact form for the translation matrix K^'n'^i dqpends on the choice of radial basis functions used in Eqn. 9. Eqn. 19 has been used previously to efficiently derive a new set of translated expansion coefficients from a reference set of expansion coefQcients as described in Danos, M., and Maximon, L. C, "Multipole matrix elements of the translation operator", J. Math. Phys.,6(l),766-778,1965; Tahnan, J. D., "Special Functions: A Group Theoretical Approach", W. A. Benjamin Inc., New York, 1968; all of which are hereby incorporated by reference in their entirety. In one embodiment, the entire set of translation matrix elements, K^',j'|„|, niay be pre-coraputed for all relevant values of n, n', 1,1', & m for a finite order of expansion, N, and stored on a computer readable and recordable medium for future retrieval as needed. This is ivantageous since calculation of the overlap integrals which define each translation naatrix lement can be very costly and yet for a given finite order of expansion and a given form for le radial basis vjs used in Eqn. 8, the calculations need to be done only once and the isultant K-matrix is applicable to any molecular subset regardless of size and shape. Moreover, for the large values of N, the number of translation matrix elements is 0(N'). In Fig. 3, in step 320, the method continues with providing a first rotation opwator Bpresenting rotation (change of orientation) of the coordinate based representation 810 of the atCTnal and external volume functions for the first molecular subset 410 with respect to the iJartesian frame 830 co-located with the molecular center 850 of molecular subset 410 in the oint coordinate system. The term "rotation operator" generally refers to an operator that, when applied to a )oint, results in the point's rotation about an axis as defined by the rotation oporator. The )perator can be applied to any collections of points, e.g., a line, a curve, a surface, or a /olume. As described with regards to Eqn. 3, any rotation in 3-D can be represented by a set Df three Euler angles. In one embodiment, different orientations of the coordinate based representation 810 of the internal and external volume functions for the first molecular subset 410 wifli respect to the Cartesian frame 830 are generally represented by a set of three Eulor angles representing roll (ai), pitch (Pi), and yaw (yO, as shown in Fig. 8. In another wnbodimcnt the roll angle (aj), describing rotation with respect to the z-axis of the Cartesian firame 830, is ignored since with respect to the common z-axis of the joint coordinate system only the relative orientation between the two molecular subsets is relevant. Then the orientation of molecular subset 410 with respect to Cartesian frame 830 is fully described by a pair of Buler angles (Pi, jj). In another embodiment the angles need not be Euler angles and in fact depend on the choice of joint coordinate systan. In one embodiment, the first rotation operator is a matrix function of (ai,pi, YI). Then the first rotation operator can be directly applied to the set of reference expansion coefficients for the intemal and external volume functions of molecular subset 410. In one embodimwit, using the joint (R,pi,yi,a2,P2,Y2) coordinate system of Fig. 8, in conjunction with the radial/spherical harmonics expansion of Eqn. 9, the first rotation operator representing the rotation of the coordinate based representation of the intemal and external volume functions for molecular subset 410 from (a]=0,Pi=0, yi=0) to arbitrary (aj.Pi, yO is directly applied to le reference set of expansion coefficients )^l^ } for the internal and external volume unctions for molecular subset 410 according to the following rule: vho-e (ajf" I are the new rotated set of expansion coefficients, (n, I, m) are the quantum nmibers for the new rotated »cpansion coefficient, m' denotes the magnetic moment quantum lumber for one of the set of reference expansion coefficients, )^^J^ j, the summation is over ill possible values of m', and R'^n* ^® matrix elements of a block diagonal matrix such that sach R^'Uenotes a (21+1)*(21+1) block sub matrix. This property tiiat the harmonic t i sxpansion coefficients transform amongst themselves under rotation in a similar way in which rotations transform the (x, y, z) coordinates in Cartesian frame was first presented in the context of molecular shapes by Leicester, S. E., Firmey, J. L., and Bywatar, R. P., in "A Quantitative Representation of Molecular-Surface Shape. 1. "nieory and Development of the Method", (1994), J. Mathematical Cbemistr)', 16(3-4), 315-341; all of which is hereby incorporated by reference in its entirety. For an arbitrary Euler rotation with angles (a,p,Y) and for a pair of positive magnetic moment quantum numbers, m and m', the individual matrix elements are computable in tcmis of Wigner rotation matrix elements, d mm'(P). as follows: where d mm' (P). the elements of the Wigno- rotation matrix are giv«i by: with kj = max (0, m-m'), kj = minO - m', 1 + m), and C(l,m,k) being a constant function. Similar forms exist for the other eight possible signed pairs of m and m'. For further details on Wigner matrix elements, refer to Su, Z., and Coppens, P., J. Applied Crystallography, 27, 89-91(1994); all of which is hereby incorporated by reference in their entirety. In one embodiment, where (aj=0) for all rotations of molecular subset 410, Eqn. 21 simplifies and the R'HOT' matrix elements are functions of (P i, yi) alone. For basis expansions other than the radial / spherical harmonics expansion of Eqn.9, ins. 20 and 21 will be replaced by appropriate analogs depending on the choice of angular isis functions, with suitable indices representing each basis function. In Fig. 3, also in step 320, the method continues with providing a second rotation perator representing rotation of the coordinate based represcaitation 820 of the internal and xtemal volume functions for the second molecular subset 420 with respect to the Cartesian •ame 840 co-located with the molecular center 860 of molecular subset 420 in the joint oordinate system. As with the first molecular subset 410, difforent orientations of the oordinate based representation 820 of the internal and external volume functions for the econd molecular subset 420 with respect to the Cartesian frame 840 are generally epreseated by a set of three Euler angles representing roll (aj), pitch (PJ). and yaw (ya), as ihown in Fig. 8. In another embodiment the angles need not be Euler angles and in fact lepend on the choice of joint coordinate system. In one embodiment, the second rotation operator is a matrix function of (a2, Pa, 72). rhen the second rotation operator can be directly applied to the set of reference expansion coefficients for the internal and external volume functions of molecular subset 420 in a marmer similar to the application of the first rotation operator to the set of referoice expansion coefficients for the internal and external volume functions of molecular subset 410. In one embodiment, the matrix function representing the second rotation operator can be split up into two distinct rotation operators, the first being a function of (P2,72) alone (i.e., a2=0) and the second being a function of the roll Euler angle, aa, alone (i.e., (Pa^O, 72=0)). Thus either of these two rotation operators can be applied first to the reference set of expansion coefficients in order to obtain an intermediate rotated set of coefficients and the remaining operator then applied in succession in order to generate a final resultant set of rotated coefficients. In such an embodiment, the two rotation operators are designated as the second and third rotation operators in order to avoid confusion regarding the first rotation i operator for molecular subset 410. Moreover, in this embodiment, when in conjunction with j the radial / spherical harmonics expansion of eqns. 9,20, and 21 can be applied for determining the result of application of each rotation opaator to the second molecular subset 420, in which case the application of the third rotation operator reduces to simple multiplication by constants and sines and cosines of the quantity (m'a). In amother embodiment, similar to the work of Ritchie et al, the simplified form for le third rotation operator permits direct application of the third rotation operator to omputed shape complementarity scores themselves, as described below, as opposed to atermediate rotated expansion coefficients for the volume functions associated with the econd molecular subset 420. In Fig. 3, in step 3 22, after the translation operators are defined, sets of translated ixpansion coefficients are constructed for the first molecular subset 410 from the sets of eference expansion coefficirents for the internal and external volume functions of molecular subset 410. The term "translated expansion coefficiattts" generally refers to a set of jxpansion coeffici«its obtained by applying a translation operator to another set of expansion :5oefBcients. As discussed above, step 310 provides for an axial sampling scheme Comprised of axial sample points which delimit the allowed values of the intermolecular separation, R, in Fig. 8 as applied to the relative translation of the two molecular subsets. In order to account for all allowed relative translations of the two molecular subsets, it is necessary to compute a set of translated expansion coefficients for both the internal and external volume functions of the first molecular subset 410, |aj^ (R = R,)), corresponding to each distinct axial sample point, Rj, in the axial sampling scheme. As discussed above, this is accomplished via direct application of a translation operator in the form of a matrix multiplication to the referaice sets of expansion coefficients for the first molecular subset 410, [aJ^{R = 0)}. In one embodiment, where the radial / sphttical harmonics expansion of Eqn. 9 is utilized, Eqn. 19 governs the construction of rnim(f^ = Rj) j for all axial sample points. Any and all permutations of the order in which axial sample points are visited is permitted, so long as in the «id the construction is completed for all axial sample points. In another embodiment, molecular subset 410 is held fixed, and the translation operator is directly applied instead to the reference sets of expansion coefficiwits for the internal and external volume functions of the second molecular subset 420, ^1^ (R = 0)}. Since only relative translation of the two molecular subsets in meaningful, it is necessary to apply the translation operator to the coordinate based representations for T and a for only one of the two molecular subsets. In Fig. 3, in step 324, after the rotation operators are defined, sets of rotated expansion lefficients are constructed for the second molecular subset 420 from the sets of reference Lpansion coefficients for the internal and external volume functions of molecular subset 20. The term "rotated expansion coefficients" generally refers to a set of expansion ^efficients obtained by applying a rotation operator to another set of expansion coefficients. iS discussed above, step 312 provides for a second spherical sampling scheme comprised of pherical sample points which delimit the allowed values of the pitch and yaw Euler angles, h* 72), in Fig. 8 as applied to orientation of the second molecular subset 420. Also as liscussed above, step 314 provides for an angular sampling scheme comprised of angular lample points which delimit the allowed values of the roll Euler angle, ai, in Fig. 8 as ipplied to rotation of the second molecular subset 420 with respect the joint z-axis. In order to account for all allowed orientations of the second molecular subset 420, it is necessary to compute a set of rotated expansion coefficients for both the intemal and external volume functions of the second molecular subset 420, { b^^ (ai = ay, p2 = p2j, Y2= Tak)}* corresponding to each distinct angular sample point, azu in the angular sampling sch^ne and each distinct spherical sample point, (Pzj, Y2k)> in the second spherical sampling scheme, i.e., (a2i,P2j, 720 e Cartesian product of the angular sampling schemne and the second spherical sampling scheme. As discussed above, this computation is accomplished via direct application of a rotation operator in the form of a matrix multiplication to the refCTence sets of expansion coefficients for the second molecular subset 420, pJJ^ j. In one embodiment, where the radial / spherical hamionics expansion of Eqn. 9 is utilized, Eqn. 20 governs the construction of {t>ta («2 == a2i, p2 = p2j, Y2 *= y2k)} for all (a2i,p2j, Y2k) e Cartesian product of the angular sampling scheme and the second spherical sampling scheme. Any and all permutations of the order in which the orientations (a2i,P2j, 72k) are visited is permitted, so long as in the end the construction is completed for all permitted (a2j,P2j. 72k)- Also as discussed above, in one embodiment the construction can be accomplished by two distinct rotational operators, the first a function of the pitch and yaw Euler angles, (P2,72), and the second a function solely of the roll Euler angle, a2. Moreover, in another embodiment, similar to the work of Ritchie et al, the latter operator (designated previously as the "third rotation operator") can be deferred until generation of shape complementarity scores, as described below. In Fig. 3, in step 326, sets of transformed expansion coefficients are constructed for 5 first molecular subset 410 from the sets of translated expansion coefficients generated in g. 3, step 322, for the internal and external volume functions of molecular subset 410. The ■m "transformed expansion coefficients" generally refers to a set of expansion coefficients •tained by applying an operator representing an arbitrary linear transformation on another t of expansion coefficiaits. This linear transformation may be the composition of one or ore translation and / or rotation operators. As discussed above, st^ 312 provides for a first spherical sampling scheme )mprised of spherical sample points which delimit the allowed values of the pitch and yaw uler angles, (pj, 71), in Fig. 8 as applied to orientation of the first molecular subset 410. As iscussed above, each set of translated expansion coefficients corresponds to an axial sample oint of an axial sampling scheme which delimits the allowed values of the intermolecular eparation, R, in Fig. 8 as applied to the relative translation of the two molecular subsets, hi irder to account for all allowed configurations (both relative orientations and translations) of he first molecular subset 410 relative to the second, it is necessary to compute a set of ransformed expansion coefficients for both the internal and external volume functions of the 3rst molecular subset 410, {Snlm (R ~ Ri, ai = 0, Pi = Pij, yi = yjk)}, corresponding to each iistinct axial sample point, Ri, in the axial sampling scheme and each distinct sphmcal sample point, (Pij,yiic), in the first spherical sampling schane, i.e., (Rj, ai = 0, Pij, yjk) e Cartesian product of the axial sampling scheme and the first sphaical sampling stheme. As discussed above, this computation is accomplished via direct application of a first rotation operator in the form of a matrix multiplication to the translated sets of expansion coefficients of step 318 for the first molecular subset 410, {3^^ (R = Ri)} • In one embodiment, where the radial / spherical harmonics expansion of Eqn. 9 is utilized, a variant of Eqn. 20 governs the constraction of { aj^^ (R = R5, aj =! 0, pi = p,j, yj = yi^)} in terms of {3 Jim (R ~ Ri)} for all (Ri, tti = 0, pij, Yik) € Cartesian product of the axial sampling scheme and the first spherical sampling scheme. Any and all permutations of the order in which (Rj, ai = 0, Pij.Yik) are visited are permitted, so long as in the end the construction is completed for all permitted (Rj, aj = 0, Py,yiu). Due to commutativity, the transformed coefficients for the first molecular subset 410 can be generated by the application of the first rotation operator and the translation operator in any order. Operations are "commutative" if the order in which they are done does not feet the results of the operations. The first rotation operator commutes with the translation )erator, so it is possible to halve instead applied the first rotation operator to the set of ference expansion coefficients, in order to generate sets of rotated coefficients for the ttemal and external volume functions for the first molecular subset 410, in a marmer similar I step 320 for the second molecular subset 420. However, as will be discussed below in regards to the generation of shape omplementarity scores, it is far more efficient in terms of computations (and potential torage) to generate sets of translated coefficients for one axial sample point at a time and to iien subsequently apply the first rotation operator in order to generate the sets of transformed ©efficients for the first molecular subset 410. In Fig. 3, in step 328, a shape complementarity score is defined. The shape omplementarity score represents a correlation between the internal volume function of the iirst molecular subset 410 with the external volume function of the second molecular subset 420, and a correlation between the internal volume function of the second molecular subset <^20 with the external volume function ofthe first molecular subset 410. This correlation represents the shape complementarity of the first and second molecular subsets for one relative position and orientation of the coordinate based represmtations of the first and second molecular subsets in the joint coordinate system. The shape complementarity score is computed in terms of the set of transformed expansion coefficients for the first molecular subset and the set of rotated (also referred to as "transformed") expansion coefficients for the second molecular subset, corresponding to a position and orientation of the first molecular subset and to a position and orientation of the second molecular subset. In one embodiment, the shape complementarity score, S, is defmed, in a marmer similar to that used in Ritchie et al, as follows: S = IoaTbdV + JobtadV-DjTaTbdV [Eqn. 23] where (Oa, T,) respectively refer to the external and internal volume functions of molecular subset 410, (Ob, Tb) respectively refer to the external and internal volume functions of molecular subset 420, D c: 0 is a constant, and the integrals over all 3-D space. The shape complementarity score of Eqn. 23 is based on assigning a positive value to the favorable overlap of the first external volun^e function and the second internal volume function, and a positive value to the favorable overlap of the first internal volume function and the second external volume function, and a negative value to the unfavorable overlap of the first and second internal volume functions, in the spirit of the Pauli Exclusion Principle. ence, the constant D in Eqn, 23 can be considered a penalty constant that weighs the third tegral relative to the first two. Shape complementarity scores are computed for different relative positions and dentations of the first and second molecular subsets 410 and 420, shown in Fig. 13 In this tnbodiment, as shown in Fig. 13 a high score is geneerated for a configuration when the xtemal volume function of the first molecular subset (aj) overlaps with the internal volume imction of the second molecular subset (T2) and vice-versa, such as when the molecular ubsets are oriented and positioned as indicated by reference numeral 1302, thus representing 1 more optimal fit. Lower scores are generally calculated for configurations when there is ittle or no overlap, as indicated by reference numeral 1306. Lower scores, and sometimes jven negative scores, are generated for configurations wh«i the internal volume function of ihe first molecular subset (xj) significantly overlaps with the internal volume function of the second molecular subset (ta) as indicated by reference numeral 1304. In Fig. 3, in step 328, a plurality of shape complementarity scores are generated by iterating over the set of sampled configurations for the molecular combination, where the set of sampled configurations is the Cartesian product of the set of sampled poses for the first molecular subset and the set of sampled poses for the second molecular subset. The iteration ova: the set of sampled configurations for the molecular combination, for the purpose of gaierating the plurality of shape complementarity scores, can be performed in any order. High correlations are achieved when the internal volume of the first molecular subset overlaps highly with the external volume of the second molecular subset and vice-versa. Low correlations result when there is lithe or no overlap between the volume functions from different molecular subsets or when significant overlap between two internal volume functions, representing unfavorable atomic overlaps. In one embodiment, the optimal fit (for example the predicted binding mode) is generally decided based on the particular configuration, i.e., relative position and orientation, that yields the highest shape complementarity score. In another embodiment, the magnitude of the best score, or the top x% of scores, determines the results of the analysis of the molecular combination of the two molecular subsets. In another embodiment, all shape complementarity scores below a preset numerical threshold are rejected, and only those configurations with passing scores are retained for further analysis. In yet another embodiment, the shape complementarity scores are filtered based on an adaptive threshold dependent on observed statistics of the scores as they are generated. In yet another nbodiment, the statistical analysis of both passing score magnitudes, as well as mltidimensional clustering of the relative position and orientation coordinates of passing onfigurations, is used to predict the binding mode and/or assess the nature and likelihood of le molecular combination. In yet another embodiment, the above strategies can be used when screening a oUection of second molecular subsets against the same first molecular subset 410 in order to iredict potential binding modes and estimate binding affinity based on computations of shape omplementarity, in order to select promising candidates for further downstream processing n the drug discovery pipeline. In one embodiment, the plurality of shape complemwitarity scores is calculated at one ^alue for the order of the expansion, Ni, and then the results are quantitatively analyzed iccording to certain decision criteria. In another embodiment, the decision criteria are based Dn a cluster analysis of the shape complementarity scores. As used herein, the term "cluster analysis" generally refers to a multivariate analysis technique that seeks to organize information about variables so that relatively homogeneous groups, or "clusters," can be formed. The clusters formed with this family of methods should be highly, internally homogenous (mwnbCTS from same cluster are similar to one another) and highly, extomally heterogeneous (members of each cluster are not like members of other clusters). A further plurality of shape complementarity scores may then be calculated at a higher value for the order of the expansion, N2 > Ni, based on results of the quantitative analysis. The shape complementarity scores may be computed at the higiher expansion order, N2, only at those sample points for which the corresponding shape complementarity score computed at the lower expansion order, Nj, satisfies the decision criteria imposed by the aforementioned quantitative analysis. Generally, it is inefficient to directly evaluate Eqn. 23 in its integral fomi. By first applying a basis expansion to the coordinate based representations of the internal and external volume functions of the two molecular subsets to obtain reference sets of expansion coefficients and then using appropriate translation and rotation operators to generate sets of transformed expansion coefficients for the first molecular subset 410 corresponding to the sampled poses of the first molecular subset, and to likewise generate sets of transformed expansion coefficients for the second molecular subset 420 corresponding to the sampled poses of the second molecular subset, the shape complementarity score for a given configuration of the molecular combination can be computed efficiently and to arbitrary precision based on the magnitude of N, the order of the expansion. In one embodiment, within the context of the joint coordinate system of Fig. 8 and ling the radial / spherical harmonics expansion of Eqn. 9, and eqns. 18,19, and 20 to mstruct and transform reference sets of expansions coefficients for the coordinate based jpresentations of both the internal and external volume functions for both molecular subsets 10 and 420, Eqn. 23 can be rewritten as follows: /here ^„/„, = i>nim ~ ^nmi. the transformed expansion coefficients for the first molecular ubset are evaluated at the sample point {ajjj^ (R = Rj, aj = 0, Pi = pij, YJ = yut)}, and the otated expansion coefficients for the second molecular subset are evaluated at the sample In the embodiment where the third rotation operator is directly applied to the computed shape complementarity scores themselves, the score is computed in two steps. In the first step, two intermediate factors Aj^ and A",, are computed, where m denote the negative values of m, the transformed expansion coefficients for the first molecular subset are as before {aj^^ (R = Ri, ai = 0, Pi = Pij, yj = yjj.)} but the rotated expansion coefficients for the second molecular subset are at the sample point, { b JjJ^ (aj = 0, p2 = Pajj Yz *= Yzk)} • In the second step, the score is given by, 26] where m is the azimuthal quantum number, and aa represents the third rotation operator. Splitting off the third rotation operator in the above marmer reduces the total computation significantly. As described above, a plurality of shape complementarity scores is generated for each Listing element of the set of sampled configurations of the molecular combination, and can 5 generated in any order. In Fig. 8, this represents a sampling of shape complementarity jores over a six-dimensional space representing the relative positions and orientations of the vo molecular subsets as given by {(R == Ri, Pi =» Pij.Yi = Yik> a2'= o-a, p2= hm, 72== Y2n)} ?here (Ri) refers to the elements of the axial sampling scheme of step 310, {pij.Yik} to the lements of the first spherical sampling scheme of step 312, {p2m,Y2n} to the elanents of the econd spherical sampling scheme, (also of step 312), and {aji} to the elements of the ingular sampling scheme of step 314. For reasonable sampling resolution for each sampling scheme, the total number of thape complementarity scores can be very large. For example, if there are 50 axial sample joints, each 1 A apart, 1000 first spherical sample points ftom an icosahedral mesh, 1000 second spherical sample points from an icosahedral mesh, and 100 angular sample points, this represents approximately five billion scores. However, reduction of the sampling resolution can lead to unacceptable inaccuracies m the final prediction and characterization of the optimal binding mode for the two molecular subsets. Thus efficiency in performing the repeated computation of Eqn. 24 for S=S(R, PI,YI, OLI, PLYI) for each sampled configuration is important, whether it be accomplished via means of computer software and/or hardware. In other embodiments, as discussed below, the method provides for fiather increased computational efficiency when considering the screening of a collection of second molecular subsets against the same first molecular subset 410. In order to increase efSciency when there is more prior knowledge about the first molecular subset 410, in one embodiment, the computation of shape complementarity scores is restricted to a subset of the possible orientations of the first molecular subset by constraining the first spherical sample points to a subset of the surface of the unit sphere. In another embodiment, the computation of shape complementarity scores can be restricted to a subset of the possible orientations of the first molecular subset by placing limits on the pitch and yaw Euler angles for the firet molecular subset. In another embodiment, when the first molecular subset 410 includes a biopolymer with one or more known active sites, the computation of shape complementarity scores is restricted to a subset of possible orientations of the first molecular subset by constraining the first spherical sample points to those that lie with the active site. In the context of large-scale screening, often little prior knowledge is known about the nding kinetics of the second molecular subset, however, if prior knowledge is available, in le embodiment the computation of shape complementarity scores is further restricted to a ibset of the possible orientations of the second molecular subset by placing limits on the igular sample points for the second molecular subset and/or placing limits on the roll, pitch, id yaw Euler angles for the second molecular subset In order to increase computational efficiency regardless of prior knowledge, previous 'ork, such as that of Ritchie et al, has employed a strategy as shown in Fig. 14. Step 1420 orresponds to the direct ^q)lication of the first rotation operator to the set of reference xpansion coefficients for molecular subset 410 m ordw to gen«"ate a set of rotated Ax.o oefficients, i.e., { "nlm (Pi = Pij. Yi = Yuc)} at each distinct (Pij.Yik). Optionally, step 1425, acre relevant in the context of large-scale screening, shows the complete set of rotated ioefficients for molecular subset 410 generated in 1420, corresponding to all first spherical '.ample points (Pij,Yik), bdng subsequently stored on a computer readable and recordable nedium. Step 1440 shows the application of the second rotation operator to the set of reference expansion coefficients for molecular subset 420 in order to generate a set of rotated b *T,CT w„>...w.w«.o,..«., ^ nim (p2= p2j, Yz^ Y2k)} at each distmct (p2j,Y2k). Optionally, step 1445, shows the complete set of rotated coefficients for molecular subset 420 generated in step 1440, corresponding to all second spherical sample points (p2j,Y2k)i being subsequently stored on a computer readable and recordable medium. Step 1450 then shows the entire set of shape complementarity scores being constructed for one specifically chosen axial sample point, Rjo, in the following marmer. First all of the rotated expansion coefficients for molecular subset 410 are retrieved fi-om the storage medium of step 1425 in step 1460 (optional, only needed if step 1425 was performed). Then the translation operator conesponding to a displacement by R = Rio, is ^plied to each and every retrieved set of rotated coefficients in step 1470, generating corresponding sets of transformed coefficients for the internal and external volume functions of molecular subset 410, i.e., { ^nlm (R = Rio, Pi = Pij, Yi ~ Yik)} for one axial sample point, Rio, and for all first spherical sample points (Pij.Yik). Then all of the rotated expansion coefficients for molecular subset 420 are retrieved fi-om the storage mediian of step 1445 in sp 1480 (optional, only needed if step 1445 was performed, but if not then step 1450 must ill occur before proceeding). Then continuing with the description of step 1450, the set of transformed coefficients om step 1470 and the set of rotated coefficients from step 1480 (or step 1440 if no storage) re combined in step 1485 in order to compute shape complementarity scores corresponding ) {S(ioj.k.ro^) = S(R = Rio, PI = Pij, Yi = Yik, h - hm Yz = Y2n)} for one axial sample point, Rjo, 11 first spherical sample points (p2j.Y2k)> and all second spherical sample points (P2j,Y2k). ccording to Eqn. 24 or an alternative form such as in Eqn. 25. In step 1490, only necessary f step 1485 utilized Eqn. 25 as opposed to 24, the third rotation operator is applied directly to he scores obtained in step 1485, according to Eqn. 26, in order to construct further scores »rresponding {S(joj.k,i.in.n) = S(R = Rjo, pi = Pij, YI * Yik, Oa* aii, P2 = p2m, Y2=Y2n)} for one ucial sample point, Rio, all first spherical sample points (P2j,Y2k), all second spherical sample points (P2j,Y2k), and all angular sample points (a2i) descaibing rotation of molecular subset 420 around the z-axis of the joint coordinate system. In step 1495, the resultant scores from step 1490,or step 1485 if 1490 was skipped, are delivered for thresholding and passing scores archived for future examination. The entirety of step 1450 (including steps 1460,1470,1480, 1485,1490, and 1495) is then repeated for another distinct axial sample point, Rii, and so on for all axial sample points in {Ri}. In the above procedure outlined in Fig. 14, step 1490 can be skipped if instead the third rotation operator was included in step 1440 as previously discussed. Moreover, steps 1425 and 1445 are optional, in that the coefficients need not be stored (and thus not retrieved so steps 1460 and 1480 become superfluous), in which case all subsequent steps which use the retrieved sets of coefficients of steps 1460 and 1480 must instead have the sets of coefficients ready at hand. However, foregoing stqps 1425 and / or 1445, is generally impractical since for order of expansion, N »1, there are simply to many coefficirats to allow immediate ('orrthe-fly') computation of ail the coefficients in a generic software and/or hardware pipeline. It is also possible to apply variants of the Fig. 14 steps that operate on subsets of the sets of expansion coefficients at each stage, as opposed to the complete sets themselves, but overall this will generally not significantly impact the performance of the procedure in Fig. 14. The procedure in Fig. 14 has a number of computational drawbacks. As will now be detailed, these drawbacks will also have major detrimental impact on any hardware architecture when the order of expansion, N, is large, i.e., N » 1. First, the number of pansion coefficients for one volume function for one molecular subset scales as 0(N^) in rms of the order of the expansion, N. For example for N=25, there are 5525 distinct ►efficients, in the set of reference expansion coefficients {^im) for one volume function for le molecular subset. Thus Eqn. 24, which minimally requires two sets of transformed )efficients for molecular subset 410 (one for x, and one for o,) and two similar rotated sets ►r molecular subset 420, there are, for N = 25,22,100 coefficients involved in the dculation of just one shape complementarity score. Granted, a reasonable ordering, when nposed on the iteration of shape complemmtarity scores through the six-dimensional space f Fig. 8, allows for keeping half the coefficiaits fixed {e.g., those for molecular subset 410) /hile computing or retrieving the other half (e.g., those for molecular subset 420) for each ew shape complementarity score, but this still represents a significant amount of monory or 'O bandwidth. Moreover, the application of Eqn. 20 to just one set of reference coefficients (i.e., ©rresponding to one spherical sample point) is 0(N*) operations, whereas the application of iqn. 19 to just one set of reference coefficients (i.e., corresponding to one axial sample point) s 0(N^) operations. Assuming that there are N,phere,i first spherical sample points, Ns,*efe,2 second spherical sample points, and NR axial sample points, and that the rotated coefficients for molecular subsets 410 and 420 will be stored on a computer readable and recordable nedium and later retrieved, the procedure outlined in Fig. 14 requires 2*Ng,Aere.i*0(N*) operations to gwierate the rotated coefficients for molecular subset 410 in step 1410, 2*Nsphete^*0(N^) operations to generate the rotated coefficients for molecular subset 420 in step 1440 and 2*Nspheo!,i*0(N^) operations in step 1470 to generate the transformed coefficients necessary for computing the shape complementarity scores in step 1485. The latter stqp in 1470 will be repeated NR times over the course of the entire process leachng to 2*Nspbere,i*NR*0(N*) Operations for translation of coefficients while steps 1420 and 1440 are performed only once assuming the storage and retrieval stqps of 1425,1445,1460, and 1480, bringing the total count to 2*(N,phere.i + Nsphc«a)*0(N^ + 2*N,pbe«.i*NR*0(N^). For example, forN^here.i ^Nsphere^" 1000, NR= 50, andN=25, this corresponds to approx. lO""'^ calculations. In the present example, this is comparable to the number of operations specific to the actual score computation represented by eqns. 25 and 26 in steps 1485 and 1490. For Nang angular sample points for the third rotation operator, the cost of the score calculations alone is Nsphere,i*Nsphere^*NR*[4*0(N^) + 2*N*N,ng]. This is roughly 1.3 x 10'^ operations for Nang = 72 angular sample points (i.e., Aaj « 5"). For the purposes of hardware architecture, the storage and bandwidth requirements of e procedure outlined in Fig. 14 are even more costly. In fact, the full set of all rotated (efficients for molecular subsets 410 and 420 represent a significant amount of data. For :ample for Nsphet«,i = Nsphere^ = 1000, N=25, assuming 32 bits precision for each value, this normts to more than 700 Mbits, a value generally far too large to store orrchip in registers rSRAM. Ifinstead the data is stored on DRAM or equivalent memory off-chip, the squirements for the memory bandwidth in the retrieval steps 1460 and 1480 are enormous id impractical, unless the hardware pipeline devoted to step 1450 is purposefully slowed to crawl; also impractical. The prospect of storing the rotated coefficients on an I/O device, ach as a disk, is even less appealing due to the low rates of I/O bandwidth relative to lemory bandwidth available for general computer and / or specialized hardware systems, liis also says nothing of the large amount of data n^esented by the precalculated values for he translation matrix elements, K^'„'|^, of Eqn. 23 in step 1470, even exploiting various lymmetries on the quantum numbers. If alternatively, when screening molecular subset 410 against a collection of nolecular subsets 420, all sets of transformed coefficients for molecular subset 410 were pre-X)mputed and stored off-chip, this would represent roughly 1 billion coefficient values that must be accessed from monory or an I/O device in the presoit example of Nsphere.i = 1000, NIR^ 50, and N=25, for just processing just one pair of molecular subsets; very impractical. The above drawbacks in terms of both operational cost and especially the storage and retrieval requirements are generally restrictive for any hardware architecture and will in fact severely limit the efficiency of any general computer software and/or hardware system. In order to increase computational efficiency regardless of prior knowledge and to eliminate the drawbacks discussed above, one embodiment of the current invention employs a strategy as shown in Fig. 15. Step 1520 corresponds to the direct application of the translation operator to the set of reference expansion coefficients for molecular subset 410 in order to generate a set of translated coefficients, i.e., {"nim (R = Rj)} at each distinct axial sample point, Rj. Step 1525 shows the complete set of translated coefficients for molecular subset 410 generated in 1520, corresponding to all axial sample points, (Rj), being subsequently stored on a computer readable and recordable medium. In Fig. 15, step 1530 then shows the entire set of shape complementarity scores being constructed for one specifically chosen axial sample point, Rjo, in the following maimo-. First all of the translated expansion coefficients for molecular subset 410 corresponding to the hosen axial sample point, Rjo, are retrieved from the storage medium of step 1525 in step 540. Then, in step 1550, the second rotation operator is applied to the set of reference ixpansion coefficients for molecular subset 420 in order to generate a set of rotated ;oefficients, i.e., {b' ^'''nim(p2 - Pzjo, Yz = Y2ko)) at a given first spherical sample point distinct P2jo,')'2ko). Then, in step 1560, the first rotation operator corresponding to a specific change in mentation by (Pyo, Yiko), is applied to each set of translated coefficients retrieved in st^ 1540, generating corresponding sets of transformed coefficients for the internal and external k'olume functions of molecular subset 410, i.e., { "nim (R = Rjo. Pi - Pijo. Yi - Yiko)} for one axial sample point, Rio, and for one specific first spherical sample point (Pijo, Yiko). Then continuing with the description of step 1530, a set of transformed coefficients from step 1560 and a set of rotated coefficients from step 1550 are combined in step 1570 in order to compute a single shape complementarity score corresponding to {S(io jo.ko,mo.no) *= S(R = Rio, Pi = Pijo, Yi = Yiko, p2= p2ino, Y2 = Yzno)} for one axial sample point, Rjo, one first spherical sample point (P2jo,Y2ko), and one second spherical sample point (PijcYaw)), according to Eqn. 24 or an alternative form such as in Eqn. 25. In step 1580, only necessary if step 1570 utilized Eqn. 25 as opposed to 24, the third rotation operator is applied directly to the scores obtained in step 1570, according to Eqn. 26, in order to construct a set of scores corresponding to {S(iojo,ko,i.niOflO) = S(R = Rw, Pi = Pijo, Yi" Yiko, aj » aji, P2 = P2mo, Y2 = Y2no)} for one axial sample point, Rjo, one first spherical sample point (P2jo,Y2ko), one second spherical sample point (p2jo,Y2koX and all angular sample points (aai) describing rotation of molecular subset 420 around the z-axis of the joint coordinate system. In step 1590, in immediate succession, the resultant scores from step 1580, or step 1570 if 1580 was skipped, are delivered for application of various decision criteria as described above. Steps 1550,1560,1570,1580, and 1590 are then repeated multiple times in order to generate the entire set of shape complementarity scores for one specifically chosen axial saiiiple point, Rjo, i.e., {S(ioj,k4.m^) = S(R = Rw, pi = Pij, Yi == Yifo ^2 = a2i, p2 = P2in, Y2 = Y2n)} for one axial sample point, R50, all first spherical sample points (P2j,Y2k), all second spherical sample points (P2j,Y2k), and all angular sample points (a2i). The number of repetitions of steps 1550 and 1560 depends on the order in which the individual scores in step 1570 are computed. The entirety of step 1530 (including steps 1540,1550,1560, 1570,1580, and 1590) is then repeated for another distinct axial sample point, Rn, and so on for all axial sample points in {Rj}. In one embodiment, step 15S0 can be skipped if instead the third rotation operator IS included in step 1550 as part of a composite rotation matrix, as previously discussed, oreover, in another embodiment, steps 1550 and 1560 may instead calculate sets of efficients for more than one spherical sample point at a time, depending on the usage of imputer readable memory to store intermediate results. In another embodiment, since it is impractical to perform steps 1550 and 1560 one me for each score generated in 1570, steps 1550 and 1560 are performed concurrently A mes and pipelined in front of step 1570, in order to feed the input requirements for enerating A* A scores in step 1570. In another embodiment step 1550 is performed A times nd step 1560 is performed B times, the results of which are stored in an intermediate omputer readable memory and pipelined in front of step 1570, in order to feed the input equirements for generating A*B scores in step 1570. In yet another embodiment, step 1540 is performed in a pipelined fashion using an ntermediate computer addressable memory so that the sets of translation coefficients X)rresponding to the next axial sample point are read in concurrently while performing one pass of step 1530. In another embodiment, for the purposes of hardware architecture, the aitire set of translated coefficients generated in step 1520 are directly stored in orrchip computer readable manory in step 1525 as opposed to off-chip computer readable mranory. In another embodiment, step 1530 performs the calculation of shape complementarity scores for more than one axial sample point in parallel and in a concurrent fashion. In anotho" embodiment step 1525 is skipped, and the translation operator is directly applied before the initiation of one pass through step 1530. The benefits of the procedure in Fig. 15 manifest thanselves in tarms of both reduced number of operations and reduced storage and memory or I/O bandwidth. In comparison to the procedure outlined in Fig. 14, assuming that there are N^ere,i first spherical sample points, Nsphere,2 secoud spherical sample points, and NR axial sraiple points, and that &e translated coefficients for molecular subsets 410 and 420 will be stored on a computer readable and recordable medium and later retrieved, now 2*NR*0(N^) operations are needed to perform step 1510 and 2*(Nspb«,.i + Nsphere^)*0(N*) operations to perform steps 1550 and 1560 for one pass through step 1530. This brings the total count to 2*NR*(Nsphere,i + Nsphere^)*0(N^) + 2*NR*0(N^), wUch Will be much smaller for large values of the expansion order N. For example, for Nsphere,! "Nsphere.2 -1000, NR= 50, and N=25, this corresponds to less than 10'^ calculations, in comparison to the lO' '''^ of the procedure in Fig. 14. In )ftware, this translates into reduced run time, and in hardware translates into either reduced irrtime and/or reduced the area, which can affect chip yields and power requirements. But far more striking is the reduction in terms of both memory storage and memory nd I/O bandwidth, especially critical for hardware designs. The entire set of translated oefficients for all axial sample points in the present example with NR = 50, N»25, and 32 bit iredsion is < 18 Mbits as opposed to the 700+ Mbits used up by the rotated coefficients of he procedure in Fig. 14. As discussed above, this data representing the entire set of ranslated coefficients can either be stored off-chip or even stored in an orrchip memory for aster access. This amount of data represents a very small comparative memory bandwidth, 5ven if stored off-chip. Also the memory or I/O bandwidth requirements for retrieving pre-jomputed translation matrix elaments are significantly reduced. Moreover, the effective use jf pre-computation of coefficients, especially valuable when screening one molecular subset against a series of other molecular subsets, becomes viable. Lastly, efficient, pipeline-able hardware architectures become possible, since the required amount of storage in orrchip memory and the amount of memory bandwidth are within the capabilities of current memory and bus technology. Embodiments of the present invention provide for pre-computation of translation expansion coefficients. These translation expansion coefficients can then be stored on a computer-readable medium. Then before computing shape complementarity scores, the stored translated expansion coefficients can be retrieved as needed for each distinct axial sample point, corresponding to a different relative translation of the two molecular subsets. In the context of assessing a likelihood of molecular combination for a pair of molecular subsets, the translation can be applied once to a set of reference coefficients for molecular subset 410 for a finite number of translational values, performed off-chip, and the results stored for subsequent use in screening against a series of second molecular subsets 420 selected from a molecule library or other collection. The molecule library is generally a database, plurality of databases, or other storage media in which a plurality of digital representations of molecular subsets are stored. Those skilled in the art should be aware that the methodology described above is applicable to a wide variety of correlatiorrbased score calculations. In addition to the volume-based shape complementarity score calculations described above, the methods are equally applicable to many other correlatiorrbased score calculations based on volumetric functions associated with molecular subsets in the context of a higih density search over relative orientations and positions of the two molecular subsets. It will be understood that the above described arrangements of apparatus and the lethod there from are merely illustrative of applications of the principles of this invention id many other embodiments and modifications may be made without departing from the jirit and scope of the invention as defined in the claims. We Claim: 1. A method of determining whether a first molecular subset is a lead candidate for a target biomolecule, wherein the target biomolecule is represented as a second molecular subset, the method comprising: providing a basis expansion of volume functions of the molecular subsets by: defining a molecular surface for each molecular subset based on locations of a plurality of surface atoms of the corresponding molecular subset; generating separately first and second internal volume functions as representations of subsets of corresponding volumes enclosed by the first and second molecular surfaces; generating separately first and second external volume functions as representations of subsets of corresponding volumes external to the first and second molecular surfaces; representing each of the internal and external volume functions by a respective reference set of expansion coefficients of a basis expansion, where the first reference set is for the first internal and external volume functions and the second reference set is for the second internal and external volume functions; defining coordinate based representations for each molecular subset using separate corresponding first and second coordinate systems; placing the coordinate based representations of the first and second molecular subsets in a joint coordinate system with separate frames for each molecular subset centered at respective molecular centers of each molecular subset, with an intermolecular axis defined there between ; defining rigid body transformations and sampling schemes by: providing a translation operator, associated with an axial sampling scheme comprising a plurality of axial sample points distributed along the intermolecular axis, reflecting discretized relative translation of the first molecular subset with respect to the second molecular subset in the joint coordinate system; providing a first rotation operator, associated with a first spherical sampling scheme comprising a plurality of spherical sample points distributed on the surface of a sphere centered on the first molecular subset, providing a second rotation operator, associated with a second spherical sampling scheme comprising a plurality of spherical sample points distributed on the surface of a sphere centered on the second molecular subset, providing a third rotation operator, associated with an angular sampling scheme comprising a plurality of angular sample points distributed on a circumference of a circle orthogonal to the intermolecular axis; defining a shape complementarity score for a molecular configuration, the shape complementarity score representing a correlation between the internal volume function of the first molecular subset with the external volume function of the second molecular subset and the internal volume function of the second molecular subset with the external volume function of the first molecular subset at a given relative position and orientation of the first and second molecular subsets in the joint coordinate system; for every axial sample point from the axial sampling scheme: constructing, with a configuration data transformation engine of a computational modeling system, a set of translated expansion coefficients for the first molecular subset, corresponding to a given axial sample point, by applying a corresponding translation operator to the first reference set of expansion coefficients for both the internal and external volume functions for the first molecular subset; storing all the sets of translated expansion coefficients for the first molecular subset on a recordable media; for each axial sample point from the axial sampling scheme and with a shape complementarity engine of the computational modeling system, retrieving, from the recordable media, the corresponding set of translated expansion coefficients for the chosen axial sample point; b) constructing a set of transformed expansion coefficients for the first molecular subset, corresponding to a set of ail configurations prescribed by the Cartesian product of the chosen axial sample point and the spherical sample points from the first spherical sampling scheme, by applying the first rotation operator to the set of translated expansions coefficients; c) constructing a set of transformed expansion coefficients for the second molecular subset, corresponding to a set of all configurations prescribed by the Cartesian product of the spherical sample points from the second spherical sampling scheme and the angular sample points from the angular sampling scheme, by applying in succession the second and third rotation operators to the second reference set of expansion coefficients for both the internal and external volume functions for the second molecular subset; and d) computing the defined shape complementarity score in terms of the set of transformed expansion coefficients for the first molecular subset, and the set of transformed expansion coefficients for the second molecular subset corresponding to each of a set of sampled configurations associated with the chosen axial sample point; and wherein at least one of the shape complementarity scores is used to determine whether the first molecular subset is a lead candidate for the target biomolecule. 2. The method as claimed in claim 1, wherein the joint coordinate system used in moving or rotating the coordinate based representations of the first molecular subset and the second molecular subset relative to one another is different from the joint coordinate system used in the final calculations of the shape complementarity scores. 3. The method as claimed in claim 1, wherein: the z-axes of the separate frames for each molecular subset are aligned with the intermolecular axis, the roll Euler angle for the first molecular subset, representing a change in orientation of the first molecular subset with respect to the z-axis of the first molecular subset's frame, is disregarded, different orientations of the first molecular subset are represented by a pair of Euler angles, namely the pitch and yaw Euler angles, and different orientations for the second molecular subset are represented by a full set of three Euler angles, including the roll Euler angle for the second molecular subset, representing a change in orientation of the second molecular subset with respect to the z-axis of the frame of the second molecular subset. 4. The method as claimed in claim 1, wherein the two successive applications of the second and third rotation operators to the coordinate based representation of the second molecular subset are combined into the application of a single combined rotation operator. 5. The method as claimed in claim 1, wherein the sets of transformed expansion coefficients of the second molecular subset are generated by applying only the second rotation operator, and the resultant shape complementarity scores are subjected to a matrix multiplication representing the third rotation operator in order to generate a plurality of shape complementarity scores that represent the set of angular sample points for the second molecular subset. 6. The method as claimed in claim 1, wherein the coordinate based representation of the internal volume function of the first molecular subset is generated for a specific coordinate system, stored on a recordable medium as a set of discrete values, each discrete value representing a portion of the information representing the coordinate based representation of the internal volume function of the first molecular subset, then the stored discrete values retrieved as needed when constructing a set of reference expansion coefficients for the internal volume function of the first molecular subset, having first converted each discrete value into another value representing the corresponding portion of the information representmg the coordinate based representation of the internal volume function of the first molecular subset in a coordinate system used to define the basis expansion, the conversion being accomplished by a suitable coordinate transformation. 7. The method as claimed in claim 1, wherein the coordinate based representation of the external volume function of the first molecular subset is generated for a specific coordinate system, stored on a recordable medium as a set of discrete values, each discrete value representing a portion of the information representing the coordinate based representation of the external volume function of the first molecular subset, then the stored discrete values retrieved as needed when constructing a set of reference expansion coefficients for the external volume function of the first molecular subset, having first converted each discrete value into another value representing the corresponding portion of the information representing the coordinate based representation of the external volume function of the first molecular subset in a coordinate system used to define the basis expansion, the conversion accomplished by a suitable coordinate transformation. 8. The method as claimed in claim 1, wherein the negative correlation representing the overlap of the first and second internal volume functions in the definition of a shape complementarity score, is multiplied by a real-valued constant so that its contribution to the shape complementarity score is opposite in sign to the positive correlation representing the net overlap of the internal volume function of one molecular subset and the external volume function of the other molecular subset. 9. The method as claimed in claim 1, wherein the plurality of sets of shape complementarity scores are calculated at one value for the order of the expansion, the shape complementarity scores are quantitatively analyzed, and a further plurality of shape complementarity scores are calculated at a higher value for the order of the expansion based on results of an intervening analysis. 10. A computational modeling system for determining whether a first molecular subset is a lead candidate for a target biomolecule, wherein the target biomolecule is represented as a second molecular subset, wherein the molecular subsets are represented with a basis expansion of volume functions provided by: defining a molecular surface for each molecular subset based on locations of a plurality of surface atoms of the corresponding molecular subset; generating separately first and second internal volume functions as representations of subsets of corresponding volumes enclosed by the first and second molecular surfaces; generating separately first and second external volume functions as representations of subsets of corresponding volumes external to the first and second molecular surfaces; representmg each of the internal and external volume functions by a respective reference set of expansion coefficients of a basis expansion, where the first reference set is for the first internal and external volume functions and the second reference set is for the second internal and external volume functions, wherein a coordinate based representations for each molecular subset uses separate corresponding first and second coordinate systems and a joint coordinate system with separate frames for each molecular subset centered at respective molecular centers of each molecular subset, with an intermolecular axis defined there between; wherein rigid body transformations and sampling schemes include: a translation operator associated with an axial sampling scheme comprising a plurality of axial sample points distributed along the intermolecular axis, reflecting discretized relative translation of the first molecular subset with respect to the second molecular subset in the joint coordinate system; a first rotation operator associated with a first spherical sampling scheme comprising a plurality of spherical sample points distributed on the surface of a sphere centered on the first molecular subset; a second rotation operator associated with a second spherical sampling scheme comprising a plurality of spherical sample points distributed on the surface of a sphere centered on the second molecular subset; and a third rotation operator associated with an angular sampling scheme comprising a plurality of angular sample points distributed on a circumference of a circle orthogonal to the intermolecular axis, the system comprising: a configuration data transformation circuit configured to construct, for every axial sample point from the axial sampling scheme, a set of translated expansion coefficients for the first molecular subset, corresponding to a given axial sample point, by applying a corresponding translation operator to the first reference set of expansion coefficients for both the internal and external volume functions for the first molecular subset; a recordable media that stores all the sets of translated expansion coefficients for the first molecular subset; a shape complementarity circuit configured to, for each axial sample point from the axial sampling scheme: a) retrieving, from the recordable media, the corresponding set of translated expansion coefficients for the chosen axial sample point; b) constructing a set of transformed expansion coefficients for the first molecular subset, corresponding to a set of all configurations prescribed by the Cartesian product of the chosen axial sample point and the spherical sample points from the first spherical sampling scheme, by applying the first rotation operator to the set of translated expansions coefficients; c) constructing a set of transformed expansion coefficients for the second molecular subset, corresponding to a set of all configurations prescribed by the Cartesian product of the spherical sample points from the second spherical sampling scheme and the angular sample points from the angular sampling scheme, by applying in succession the second and third rotation operators to the second reference set of expansion coefficients for both the internal and external volume functions for the second molecular subset; and d) computing the defined shape complementarity score in terms of the set of transformed expansion coefficients for the first molecular subset, and the set of transformed expansion coefficients for the second molecular subset corresponding to each of a set of sampled configurations associated with the chosen axial sample point, wherein the shape complementarity score represents a correlation between the internal volume function of the first molecular subset with the external, volume function of the second molecular subset and the internal volume function of the second molecular subset with the external volume function of the first molecular subset at a given relative position and orientation of the first and second molecular subsets in the joint coordinate system; and a combination post-processor circuit that determines whether the first molecular subset is a lead candidate based on the shape complementarity scores.