In the quantum mechanical picture of matter, the interaction between molecules arises from the interaction of electrons and nuclei in one molecule with the electrons and nuclei in another molecule. In the case of small molecules, accurate quantum chemical calculations readily yield the energy of isolated molecules, and the total energy of the complex. From these energies, the interaction energy in the gas phase can be obtained as
Eint = EAB - EA - EB
In the molecular mechanical picture of matter, the interaction between molecules arises from the interaction of charge distributions and deformation of molecules. In most programs, the charge distribution in the molecule is represented by partial point charges and van der Waals spheres located in the position of nuclei. The internal deformation of the molecule is descruibed as displacement of bond lengths and angles from predefined "equilibrium" values. Despite a fundamentally different description of molecular structure and energetics, the interaction energy in molecular mechanics can be obtained similarly by subtracting the energies of unbound molecules from the energy of the complex.
Calculation of interaction strength between molecules of medium size (i.e. typical drugs) is a challenging problem for quantum mechanics. The difficulties arises from the fact that the description of the London dispersion interaction requires large basis sets, and Hamiltonians that correctly capture the electron correlation effects. The London dispersion is related to molecular polarizability, and diffuse functions (i.e. the "+" in the 6-31+G(d,p) or "aug" in the aug-cc-pVTZ basis) are generally needed to describe polarizability correctly. However, using such diffuse Gaussian functions increases the computational time, and may lead to the convergence problems in the self-consistent field stage.
Choice of a suitable quantum mechanical Hamiltonian is critical for meaningful results. Some rather advanced methods, such as CISD, are fundamentally unsuitable for the description of molecular interactions because they are not size consistent. Almost all commonly used model chemistries, such as Hartree-Fock (HF), Density Functional Theory (DFT), Møller-Plesset (MP) perturbation theory, and Coupled Clusters (CC) method describe the first order electrostatic interaction between charge distributions correctly. Advanced methods that capture the effects of dynamic electron correlation (MP2, CCSD) allow the description of London dispersion interaction. However, London dispersion is completely neglected in some computationally efficient models, including the popular Hartree-Fock theory. London dispersion is also missing in many older Density Functional Theory functionals (e.g. in the popular B3LYP) but is included as a an empirical correction in some newer functionals, such as B2PLYP-D. As a result, HF or B3LYP methods would incorrectly predict that naphthalene dimers are most stable in the T-shaped geometry while more accurate MP2 or B2PLYP-D methods predict that the favored mode of interaction involves π-π stacking in the slipped-parallel geometry. Semiempirical models do a poor job predicting the lowest-energy geometry of naphthalene dimer: AM1, PM3, and PM6 predict different T-shaped forms as global minima. The image below shows the lowest-energy slipped parallel geometry for naphthalene dimer, followed by AM1, PM3, and PM6 predictions.
Accurate description of London dispersion, as well as of second order induction-polarization energy usually requires that a large, flexible basis set of Gaussian functions be used to describe molecular orbitals. Too few Gaussians will not allow sufficient polarization, and thus the induction energy is underestimated. Minimally, one set of polarization functions, and one set of diffuse functions are recommended for the calculation of intermolecular interactions. Thus, to accurately model the induced polarization component in hydrogen bonding, one would need a Hartree-Fock calculation with the 6-31++G(d,p) basis where the two "+" signs specify that diffuse functions are used both on the heavy atoms (C, N, O, ...) and on the hydrogen. The "(d)" indicates that a d-type polarization function is placed on every heavy atom, and the "(p)" indicates that the electron density at hydrogen is allowed to be non-spherical thanks to one p-type function. However, to accurately model both the induced polarization and London dispersion in hydrogen bonding, a calculation with basis as large as aug-cc-pVQZ is preferred because this basis adds multiple polarization and diffuse functions which are needed for the description of both the polarization of charge distribution and dynamic electron correlation.
Accurate calculation of interaction energy via quantum mechanics faces another hurdle. Molecular orbitals of each monomer are described by a set of Gaussian functions. However, in practical calculations, only a small number of Gaussian functions are used such that each monomer is not perfectly satisfied with its Gaussian functions. If the monomers are far apart, this less-than satisfactory basis set coverage leads to a slightly erroneous energy for a monomer. This is normally not a big deal: almost all quantum chemistry calculations suffer from such basis set incompleteness errors, but these errors cancel largely out if one compares to the energies of two isomers of the same molecule, or energies of reactants and products. However, when the two monomers approach each other, they start borrowing some Gaussian functions from each other to better describe each other's molecular orbitals. As each monomer is now using its own, as well as its partner's basis functions, the energy of the monomer (even in the absence of any electrostatic or London dispersion interaction) is lower than that of the isolated monomer. Thus, the energy of the complex becomes contaminated with the so called Basis Set Superposition Error (BSSE). The magnitude of this error is large with the commonly used basis sets. It can be removed by performing additional calculations (Counterpoise keyword in Gaussian03), but this adds to the overall effort and time to describe molecular interactions. The magnitude of this error diminishes as the basis set is increased but the elimination of BSSE would require extremely large basis sets.
In case of protein-ligand interaction, the system is too large for the meaningful use of quantum mechanics, and variety of simplified schemes, such as molecular mechanics is used to obtain the interaction energy. In the latter case, it is helpful to think of the physical origin of the interaction because this allows to derive compact expressions for the calculation of interaction energy as a function of intermolecular distance and geometry. For example, typical force fields represent molecular electron density as a collection of atom-centered point charges and the interaction between molecules is calculated from the Coulomb law. The Coulombic interaction energy is supplemented with the Lennard-Jones potential in order to describe the Pauli repulsion between electron clouds, and London dispersion. Non-polarizable force fields do not allow for redistribution of charges in the molecule as it is approached by an ion; to compensate for the lack of induction energy, the point charges are usually enhanced by 10-20% in such force fields. Polarizable force fields allow for redistribution of charges within the molecule but are significantly slower than non-polarizable force fields. In general, most molecular mechanics force fields based on fixed point charge model describe well hydrogen bonding but fail to treat more subtle interactions, such as cation-pi, quadrupole-quadrupole, or charge-transfer interactions.