*ab initio*quality force fields from condensed phase quantum-mechanics/molecular-mechanics calculations through the adaptive force matching method

^{1}, Yang Song

^{1}and Feng Wang

^{1,a)}

### Abstract

A new method called adaptive force matching (AFM) has been developed that is capable of producing high quality force fields for condensed phase simulations. This procedure involves the parametrization of force fields to reproduce *ab initio* forces obtained from condensed phase quantum-mechanics/molecular-mechanics (QM/MM) calculations. During the procedure, the MM part of the QM/MM is iteratively improved so as to approach *ab initio* quality. In this work, the AFM method has been tested to parametrize force fields for liquid water so that the resulting force fields reproduce forces calculated using the *ab initio* MP2 and the Kohn–Sham density functional theory with the Becke–Lee–Yang–Parr (BLYP) and Becke three-parameter LYP (B3LYP) exchange correlation functionals. The AFM force fields generated in this work are very simple to evaluate and are supported by most molecular dynamics (MD) codes. At the same time, the quality of the forces predicted by the AFM force fields rivals that of very expensive *ab initio* calculations and are found to successfully reproduce many experimental properties. The site-site radial distribution functions (RDFs) obtained from MD simulations using the force field generated from the BLYP functional through AFM compare favorably with the previously published RDFs from Car–Parrinello MD simulations with the same functional. Technical aspects of AFM such as the optimal QM cluster size, optimal basis set, and optimal QM method to be used with the AFM procedure are discussed in this paper.

We thank Gerrick Lindberg and Dr. Christian J. Burnham for helpful discussions. This work was supported by ACS PRF under Grant No. 47001-AC6 and NSF CAREER Award No. CHE0748628. The computer resource for this study was provided by the National Center for Supercomputing Applications under Grant No. MRAC TG-CHE070060 and by the Boston University’s Center for Scientific Computing and Visualization.

I. INTRODUCTION

II. AFM ALGORITHM

III. TECHNICAL CONSIDERATIONS

A. Choice of functional form for the force field

B. Choice of the objective function being minimized

C. Choice of basis set

D. Choice of number of QM molecules

IV. SIMULATIONS AND RESULTS

V. CONCLUSIONS

### Key Topics

- Atomic force microscopy
- 53.0
- Ab initio calculations
- 25.0
- Density functional theory
- 12.0
- Molecular dynamics
- 12.0
- Intermolecular forces
- 9.0

## Figures

Scatter plot of forces at the B3LYP level of theory using the 6-311G (blue circles), aug-cc-pvdz (black crosses), aug(dz)-cc-pvtz (green triangles), and aug-cc-pvtz (red squares) basis sets. The axis is the B3LYP/aug-cc-pvqz result. Panel (a) shows the atomic forces, panel (b) shows net molecular forces, and panel (c) shows net molecular torques multiplied by . All the points would lie on the black diagonal lines if the different basis sets resulted in identical atomic forces.

Scatter plot of forces at the B3LYP level of theory using the 6-311G (blue circles), aug-cc-pvdz (black crosses), aug(dz)-cc-pvtz (green triangles), and aug-cc-pvtz (red squares) basis sets. The axis is the B3LYP/aug-cc-pvqz result. Panel (a) shows the atomic forces, panel (b) shows net molecular forces, and panel (c) shows net molecular torques multiplied by . All the points would lie on the black diagonal lines if the different basis sets resulted in identical atomic forces.

Scatter plot of forces at the MP2 level of theory using the 6-311G (blue circles), aug-cc-pvdz (black crosses), and aug(dz)-cc-pvtz (green triangles) basis sets. The axis is the MP2/aug-cc-pvtz result. Panel (a) shows the atomic forces, panel (b) shows net molecular forces, and panel (c) shows net molecular torques multiplied by . All the points would lie on the black diagonal lines if the different basis sets resulted in identical atomic forces.

Scatter plot of forces at the MP2 level of theory using the 6-311G (blue circles), aug-cc-pvdz (black crosses), and aug(dz)-cc-pvtz (green triangles) basis sets. The axis is the MP2/aug-cc-pvtz result. Panel (a) shows the atomic forces, panel (b) shows net molecular forces, and panel (c) shows net molecular torques multiplied by . All the points would lie on the black diagonal lines if the different basis sets resulted in identical atomic forces.

Scatter plot of forces at the BLYP/aug-cc-pvdz (red circles), PBE (green crosses), and the B3LYP/aug-cc-pvdz (blue pluses) levels of theory. The axis is the MP2/aug-cc-pvdz result. Panel (a) shows the atomic forces, panel (b) shows net molecular forces, and panel (c) shows net molecular torques multiplied by . All the points would lie on the black diagonal lines if the different methods resulted in identical atomic forces.

Scatter plot of forces at the BLYP/aug-cc-pvdz (red circles), PBE (green crosses), and the B3LYP/aug-cc-pvdz (blue pluses) levels of theory. The axis is the MP2/aug-cc-pvdz result. Panel (a) shows the atomic forces, panel (b) shows net molecular forces, and panel (c) shows net molecular torques multiplied by . All the points would lie on the black diagonal lines if the different methods resulted in identical atomic forces.

Scatter plot of forces at the MP2/aug-cc-pvdz level of theory. The axis is the QCISD/aug-cc-pvdz result. Panel (a) shows the atomic forces, panel (b) shows net molecular forces, and panel (c) shows net molecular torques multiplied by . All the points would lie on the diagonal lines if the MP2 and QCISD forces are identical.

Scatter plot of forces at the MP2/aug-cc-pvdz level of theory. The axis is the QCISD/aug-cc-pvdz result. Panel (a) shows the atomic forces, panel (b) shows net molecular forces, and panel (c) shows net molecular torques multiplied by . All the points would lie on the diagonal lines if the MP2 and QCISD forces are identical.

Scatter plots of forces on the central molecules with (blue circles), (red triangles), and (black crosses). The axis corresponds to the forces for . Panel (a) shows the atomic forces, panel (b) shows net molecular forces, and panel (c) shows net molecular torques multiplied by . All the points would lie on the black diagonal lines if the atomic forces on the central molecule were the same for each of the cluster sizes considered.

Scatter plots of forces on the central molecules with (blue circles), (red triangles), and (black crosses). The axis corresponds to the forces for . Panel (a) shows the atomic forces, panel (b) shows net molecular forces, and panel (c) shows net molecular torques multiplied by . All the points would lie on the black diagonal lines if the atomic forces on the central molecule were the same for each of the cluster sizes considered.

Scatter plot of the forces for the (red circles) and SPCFw (blue crosses). The axis corresponds to the *ab initio* MP2 forces. Panel (a) shows the atomic forces, panel (b) shows net molecular forces, and panel (c) shows net molecular torques multiplied by . All the points would lie on the green diagonal lines if the different models predict forces identical to the MP2 *ab initio* forces.

Scatter plot of the forces for the (red circles) and SPCFw (blue crosses). The axis corresponds to the *ab initio* MP2 forces. Panel (a) shows the atomic forces, panel (b) shows net molecular forces, and panel (c) shows net molecular torques multiplied by . All the points would lie on the green diagonal lines if the different models predict forces identical to the MP2 *ab initio* forces.

Site-site RDFs at the zeroth, first, and second generations obtained from AFM using MP2/aug-cc-pvdz forces. The RDFs from our final force field is also shown together with experimental results. The blue thick solid lines with error bars are experimental results from Ref. 34.

Site-site RDFs at the zeroth, first, and second generations obtained from AFM using MP2/aug-cc-pvdz forces. The RDFs from our final force field is also shown together with experimental results. The blue thick solid lines with error bars are experimental results from Ref. 34.

Comparison of the RDFs from with those from CPMD simulation of Lee and Tuckerman (Ref. 37). The experimental RDFs are also plotted.

Comparison of the RDFs from with those from CPMD simulation of Lee and Tuckerman (Ref. 37). The experimental RDFs are also plotted.

Site-site RDFs from the final force fields from the MP2, BLYP, and B3LYP methods shown together with experimental RDFs.

Site-site RDFs from the final force fields from the MP2, BLYP, and B3LYP methods shown together with experimental RDFs.

Infrared spectra from simulations with the , , and force fields shown together with the experimental results.

Infrared spectra from simulations with the , , and force fields shown together with the experimental results.

## Tables

Parameters for the force fields and for the , , and force fields. The quantities and refer to the equilibrium OH and HOH intramolecular bond length and angle, respectively. We also show the parameters for the SPCFw model (Ref. 33) for comparison.

Parameters for the force fields and for the , , and force fields. The quantities and refer to the equilibrium OH and HOH intramolecular bond length and angle, respectively. We also show the parameters for the SPCFw model (Ref. 33) for comparison.

Thermodynamic properties calculated from our final models.

Thermodynamic properties calculated from our final models.

Article metrics loading...

Full text loading...

Commenting has been disabled for this content