The atomic polarizabilities are listed for each AMOEBA atom type in Table 1 . The values are the same as those derived by Thole62 except for aromatic carbon and hydrogen atoms which have been systematically refined using a series of aromatic systems, including a small carbon nanotube (see Table 3 ). The molecular polarizabilities computed using the current model are compared to experimental values for selected compounds in Table 2 . Reducing the damping factor from Thole’s original value of 0.567 to AMOEBA’s 0.39 is critical to correctly reproducing water cluster energetics.18 On the other hand, AMOEBA’s greater damping leads to a slight systematic underestimation of molecular polarizabilities. However, given the simplicity of the model, the agreement is generally satisfactory for both average polarizabilities and their anisotropies. As described above, polarization groups are defined for purposes of treating intramolecular polarization. Typically a functional group is treated as a single polarization group. For example, methylamine is a group by itself, while ethylamine has two groups: –CH2NH2 and CH3–. The groups are specified in AMOEBA parameter files in the following format: “polarize A α 0.390 B C”, where α is the polarizability for atom type A; 0.39 is the damping coefficient in Eq (2) ; B and C are possible bonded atom types that belong to the same polarization group as atom type A.
The permanent atomic multipoles were derived for each molecule from ab initio QM calculations. Ab initio geometry optimization and a subsequent single-point energy evaluation were performed at the MP2/6-311G(1d,1p) level using Gaussian 03.78 For small molecules with less than six heavy atoms, Distributed Multipole Analysis (DMA v1.279 ) was used to compute the atomic multipole moments in the global frame using the density matrix from the QM calculation. Next, the TINKER POLEDIT program rotates the atomic multipoles into a local frame and extracts Thole-based intramolecular polarization to produce permanent atomic multipole (PAM) parameters. Thus, when the AMOEBA polarization model is applied to the permanent atomic moments, the original ab initio-derived DMA is recovered. Finally, the POTENTIAL program from the TINKER package is used to optimize the permanent atomic multipole parameters by fitting to the electrostatic potential on a grid of points outside the vdW envelope of the molecule. The reference potential for the fitting step is typically derived from a single point calculation at the MP2/aug-cc-pVTZ level. Only a partial optimization to the potential grid is used to keep the atomic moments close to their DMA-derived values while still providing an improved molecular potential. The fitting approach is also useful for molecules containing symmetry-averaged atoms of the same atomic multipole type. In this case, simple arithmetic averaging would degrade the quality of the PAM. For example, in dimethyl- or trimethylamine, all the methyl hydrogen atoms are indistinguishable and adopt the same atom type. The DMA multipole values for these atoms are somewhat different due their non-equivalence in any single conformation, and PAM derived by simple averaging would lead to a large error in the molecular dipole moment. The potential-optimized PAM, where methyl hydrogens are constrained to adopt equivalent values, will reproduce almost exactly both the ab initio potential and the molecular multipole moments. Our standard procedure is to use a molecular potential grid consisting of a 2 Å shell beginning 1 Å out from the vdW surface. The DMA monopole values are generally fixed during the potential fitting procedure.
This electrostatic parameterization protocol is particularly important for larger molecules and for molecules with high symmetry. It is known that the original DMA approach tends to give “unphysical’ multipole values for large molecules when diffuse functions are included in the basis set even though the resulting electrostatic potential is correct. A recent modification of DMA80 has been put forward to address this issue. However, in our hands, the multipoles from the modified scheme seem less transferable between conformations. The above protocol allows derivation of PAM corresponding to larger basis sets than would be practical with the original DMA method. Note this procedure is different from restrained potential fits commonly used to fit fixed atomic charge models, as the starting DMA values are already quite reasonable and the fitting can be considered as a small perturbation biased toward the larger basis set potential. The overall procedure has been extensively tested in a small molecule hydration study81 (link) and will be used in future AMOEBA parameterization efforts.
Empirical vdW parameters were determined by fitting to both gas and liquid phase properties. The gas phase properties include homodimer binding energy (BSSE corrected) and structure from ab initio calculations at the MP2/aug-cc-pVTZ level or above. Liquid properties include experimental density and heat of vaporization of neat liquids. The vdW parameters were first estimated by comparing structure and energy of the AMOEBA-optimized dimer with ab initio results, and then fine-tuned to reproduce the experimental liquid density and heat of vaporization via molecular dynamics simulation. Additional homodimers at alternative configurations, heterodimers with water, and liquid properties were computed post facto for the purpose of validation. A more generic force field atom classification for vdW parameters was enforced to ensure the transferability.Table 1 lists the common vdW atom classes used by AMOEBA, together with the corresponding vdW parameters and polarizabilities. The vdW atom classes are also used to define parameters for all of the valence potential energy terms. The parameters for bonded terms, initially transferred from MM3, are optimized to reproduce ab initio geometries and vibration frequencies. In the final parameterization step, after all other parameters are fixed, torsional parameters were obtained by fitting to ab initio conformational energy profiles at the MP2/6-311++(2d,2p) level of theory.
The permanent atomic multipoles were derived for each molecule from ab initio QM calculations. Ab initio geometry optimization and a subsequent single-point energy evaluation were performed at the MP2/6-311G(1d,1p) level using Gaussian 03.78 For small molecules with less than six heavy atoms, Distributed Multipole Analysis (DMA v1.279 ) was used to compute the atomic multipole moments in the global frame using the density matrix from the QM calculation. Next, the TINKER POLEDIT program rotates the atomic multipoles into a local frame and extracts Thole-based intramolecular polarization to produce permanent atomic multipole (PAM) parameters. Thus, when the AMOEBA polarization model is applied to the permanent atomic moments, the original ab initio-derived DMA is recovered. Finally, the POTENTIAL program from the TINKER package is used to optimize the permanent atomic multipole parameters by fitting to the electrostatic potential on a grid of points outside the vdW envelope of the molecule. The reference potential for the fitting step is typically derived from a single point calculation at the MP2/aug-cc-pVTZ level. Only a partial optimization to the potential grid is used to keep the atomic moments close to their DMA-derived values while still providing an improved molecular potential. The fitting approach is also useful for molecules containing symmetry-averaged atoms of the same atomic multipole type. In this case, simple arithmetic averaging would degrade the quality of the PAM. For example, in dimethyl- or trimethylamine, all the methyl hydrogen atoms are indistinguishable and adopt the same atom type. The DMA multipole values for these atoms are somewhat different due their non-equivalence in any single conformation, and PAM derived by simple averaging would lead to a large error in the molecular dipole moment. The potential-optimized PAM, where methyl hydrogens are constrained to adopt equivalent values, will reproduce almost exactly both the ab initio potential and the molecular multipole moments. Our standard procedure is to use a molecular potential grid consisting of a 2 Å shell beginning 1 Å out from the vdW surface. The DMA monopole values are generally fixed during the potential fitting procedure.
This electrostatic parameterization protocol is particularly important for larger molecules and for molecules with high symmetry. It is known that the original DMA approach tends to give “unphysical’ multipole values for large molecules when diffuse functions are included in the basis set even though the resulting electrostatic potential is correct. A recent modification of DMA80 has been put forward to address this issue. However, in our hands, the multipoles from the modified scheme seem less transferable between conformations. The above protocol allows derivation of PAM corresponding to larger basis sets than would be practical with the original DMA method. Note this procedure is different from restrained potential fits commonly used to fit fixed atomic charge models, as the starting DMA values are already quite reasonable and the fitting can be considered as a small perturbation biased toward the larger basis set potential. The overall procedure has been extensively tested in a small molecule hydration study81 (link) and will be used in future AMOEBA parameterization efforts.
Empirical vdW parameters were determined by fitting to both gas and liquid phase properties. The gas phase properties include homodimer binding energy (BSSE corrected) and structure from ab initio calculations at the MP2/aug-cc-pVTZ level or above. Liquid properties include experimental density and heat of vaporization of neat liquids. The vdW parameters were first estimated by comparing structure and energy of the AMOEBA-optimized dimer with ab initio results, and then fine-tuned to reproduce the experimental liquid density and heat of vaporization via molecular dynamics simulation. Additional homodimers at alternative configurations, heterodimers with water, and liquid properties were computed post facto for the purpose of validation. A more generic force field atom classification for vdW parameters was enforced to ensure the transferability.
