Skip to content

JOYCE examples: parametrization of QMD-FF

General Introduction

Welcome to the introductory JOYCE tutorials. Here, we discuss the fundamentals of JOYCE through practical examples.

Transferable vs Specific Force Fields

General-purpose force fields are designed to describe a certain group of similar molecules with the same transferable parameters. A large number of molecules is drawn and some global parameters are given that reproduce targetted experimentally-measured thermodynamic or physical properties of the molecular species. Examples of such are typically the density, vaporization enthalpy or dipole moment.

Conversely, QMD-FFs are tailored to describe a specific molecular species, ranging from organic molecules to metal-organic frameworks. Concretely, QMD-FFs are built by matching selected molecular descriptors purposely computed at QM level. JOYCE is a software that automatizes the QMD-FF parametrization procedure.

Required ingredients for a Force Field

Every FF requires defining 3 ingredients:

  • Definition of Internal Coordinates (ICs)

We must define a potential energy function that depends of the atomic positions for the MD simulation. In its most basic form, JOYCE assumes the following potential energy:

\[ U(\{r\})= E_{stretch}(\{r_\mu\}) + E_{bend}(\{\theta_\mu\}) + E_{Rtors}(\{\phi_\mu\}) + E_{Ftors}(\{\delta_\mu\}) + E_{nb}^{intra}(\{r_{ij}\}) \]

Info

It is also possible to define coupling terms involving two ICs within the JOYCE parametrization -- please refer to the Manual for further information. This is an advanced topic, so we limit the discussion to the most basic parametrizations for the tutorials.

  • Functional form of the ICs

We recommend making the following distinction between hard and soft degrees of freedom:

  1. Hard/Stiff degrees of freedom have a high frequency (> 200cm\(^{-1}\)) such that the local harmonic approximation (LHA) is enough for an accurate representation.

  2. Soft/Flexible degrees of freedom have lower frequencies should be fitted to the full rotation of the internal coordinate.

All the stiff ICs are modelled with the harmonic potential. Some dihedrals are considered soft such as a phenyl-phenyl rotation, which could span the 360\(\deg\) rotation in a MD simulation. Long range terms are advised to be discarded unless required for highly flexible molecules. In particular, Coulomb terms are always null within JOYCE. Additional LJ terms may be added avoid forbidden configurations. This will highly depend on the molecular flexibility, e.g. a benzene ring would not require these owing to its high stiffness.

In short, JOYCE assumes the following functions for the ICs:

\[ E_{stretch}(\{r_\mu\}) = \frac{1}{2}\sum_\mu^{N_{bonds}}k_\mu^s(r_\mu - r_\mu^0)^2 \]
\[ E_{bend}(\{\theta_\mu\}) = \frac{1}{2}\sum_\mu^{N_{angles}}k_\mu^b(\theta_\mu - \theta_\mu^0)^2 \]
\[ E_{Rtors}(\{\phi_\mu\}) = \frac{1}{2}\sum_\mu^{N_{Rdihedrals}}k_\mu^t(\phi_\mu - \phi_\mu^0)^2 \]
\[ E_{Ftors}(\{\delta_\mu\}) = \sum_\mu^{N_{Fdihedrals}}\sum_{j=1}^{N_{cos_\mu}}k_{j\mu}^d[1+\cos(n_{j}^\mu\delta_\mu - \gamma_{j}^\mu)] \]
\[ E_{nb}^{intra}(\{r_{ij}\}) = \sum_{i=1}^{N_{nb}}\sum_{j=i+1}^{N_{nb}} 4\epsilon^{intra}_{ij}\left[\left(\frac{\sigma^{intra}_{ij}}{r_{ij}}\right)^{12}-\left(\frac{\sigma^{intra}_{ij}}{r_{ij}}\right)^{6}\right] + \sum_{i=1}^{N_{nb}}\sum_{j=i+1}^{N_{nb}} \frac{q_iq_j}{r_{ij}},\quad q_i=q_j=0 \]
  • Force and equilibrium constants

For a given nuclear conformation, the electronic interaction energy is captured in the force constants. Then, through the potential energy function, we may solve the motion of the nuclei position though Newton's 2nd law. Hence, these force constants along with the specific functional form capture the electronic response to all the allowed nuclei displacements.

More details, including the equations employed in JOYCE, can be found in the Manual.

JOYCE Protocol

We have already discussed how the first 2 ingredients are handled by JOYCE when we introduced them. The 3rd ingredient is obtained through a Hessian minimization. Equilibrium constants are extracted from a the absolute minimum of a DFT optimization. The force constants are then fitted to minimize the difference with the QM Hessian matrix. It is possible to also fit the energy gradient, but this is not usually required. For more information regarding the methods, refer to the Manual.

The Frozen Internal Rotation Approximation (FIRA) is assumed to fit flexible internal coordinates. For that, a QM-DFT scan is required for full rotations around every the flexible ICs.

Different QM-based computational chemistry softwares, can be used to obtain the necessary data. JOYCE is compatible with both the formatted check point file format (.fchk) and the formatted file format (.fcc). QM-data obtained in other formats should be converted to compatible files.

Pyridine