Force-field model potential functions

Once the set of RIC to be employed has been chosen, the analytical form of the model potential energy functions depending on the selected RIC must be declared. In the JOYCE code, the intramolecular FF is thus expressed as a sum of different terms, namely
$\displaystyle V^{intra} = E_{stretch} + E_{bend} + E_{Rtors} + E_{Ftors} +
E_{nb}^{intra} + E_{Coupl}$     (6)

The first three terms count for the stiff IC, i.e. bond stretchings, angle bendings and stiff angle dihedrals ( $ Rdihedrals$), as those that drive the planarity of aromatic rings and are expressed with harmonic potentials:
$\displaystyle E_{stretch} =
\frac{1}{2} \sum_\mu^{N_{bonds}} k^s_\mu (r_\mu - r^0_\mu) ^2$     (7)


$\displaystyle E_{bend} =
\frac{1}{2} \sum_\mu^{N_{angles}}
k^b_\mu (\theta_\mu - \theta^0_\mu)^2$     (8)


$\displaystyle E_{Rtors} =
\frac{1}{2} \sum_\mu^{N_{Rdihedrals}}
k^t_\mu (\phi_\mu - \phi^0_\mu)^2$     (9)

Conversely, the model functions employed for soft, flexible dihedrals ( $ Fdihedrals$) are sums of periodic functions, namely
$\displaystyle E_{Ftors} = \sum_\mu^{N_{Fdihedrals}}
\sum^{N_{cos_\mu}}_{j=1}
k^{d}_{j\mu}
\big[ 1 + \cos ( n_j^\mu \delta_\mu - \gamma_j^\mu ) \big]$     (10)

where $ N_{cos_\mu}$ is the number of cosine function employed to describe the potential of the $ \delta_\mu$ dihedral. Finally, if any internal distance between atoms not directly bonded to each other is defined, the non-bonded intramolecular contribution is computed as
$\displaystyle E_{nb}^{intra}$ $\displaystyle =$ $\displaystyle \sum^{N_{nb}}_{i=1} \sum^{N_{nb}}_{j=i+1}
4 \epsilon_{ij}^{intra}...
...} } \bigg )^{12} -
\bigg(\frac{ \sigma_{ij}^{intra} } {r_{ij} } \bigg)^6 \bigg]$  
  $\displaystyle +$ $\displaystyle \sum^{N_{nb}}_{i=1} \sum^{N_{nb}}_{j=i+1}
\frac{q_iq_j}{r_{ij}}$ (11)

Note that, at difference with intermolecular energy terms reported in equations (3) and (4), $ i$ and $ j$ indexes run over atoms of the same molecule, therefore $ j$ $ \ne$ $ i$, to avoid self interaction. Furthermore, the interaction parameters $ \epsilon_{ij}^{intra}$ and $ \sigma_{ij}^{intra}$ do not have to be necessarily the same employed in the intermolecular LJ interaction reported in equation (4), although some MD software do not allow to use different LJ parameters for intra and inter-molecular interactions.

All of the aforementioned FF terms are expressed as sums of contribution, each depending on a single IC (even if, as before stated, $ r_{ij}$ can always be expressed as a function of other generalized IC's). If the last term of equation (6) is not included, such FF is often termed as diagonal. As a matter of fact the $ E_{Coupl}$ term, takes explicitly into account the coupling of two (or more) RIC. Since 2018, [21] the JOYCEcode has been equipped to handle two-variable functions designed to explicitly account for the coupling between the considered ICs. As the detailed in the original paper, [21] a more convenient formalism can be adopted to clarify the coupling functions that can be used within the JOYCEparameterization procedure. In fact, when the couling term $ E_{Coupl}$ is missing, the intramolecular energy of a purely diagonal FF can be re-written as

$\displaystyle V^{intra}(\mathbf{q}) =\sum_{a=1}^{N_{func}} p_{a} f_{a}(q_\mu)$     (12)

where $ q_\mu$ stands for any internal coordinate $ r_\mu,\theta_\mu,\phi_\mu,\delta_\mu$ reported in equations (7)-(10), $ p_a$ refers to the force constant (either $ k^s_\mu$, $ k^b_\mu$, ...), while $ f_a$ is the diagonal potential function ( $ [r_\mu - r_\mu^0]^2$, $ [\theta_\mu - \theta_\mu^0]^2$, $ [ 1 + \cos (n_j \delta_\mu - \gamma_j)]$, etc.) assigned to $ q_\mu$. Within this framework, the generalized coupling introduced in JOYCEconsists in a sum of $ N_{Coupl}$ pairwise linear terms, defined as the product between two functions, $ f$ and $ g$, each depending only on one of the coupled ICs:
$\displaystyle E_{Coupl} = \sum_a^{N_{Coupl}} p_af_a(q_\mu)g_a(q_\nu)$     (13)

Following this notation, the coupling terms between two stiff ICs, are expressed by a product of linear terms involving each considered coordinate, as for instance:

$\displaystyle E^{str-str}_{\mu\nu}(r_\nu, r_\mu) = k^{ss}_{\mu\nu}
( r_\mu - r^0_\mu) ( r_\nu - r^0_\nu)$     (14)


$\displaystyle E^{bnd-bnd}_{\mu\nu} (\theta_\mu, \theta_\nu) = k^{bb}_{\mu\nu}
( \theta_\mu - \theta^0_\mu) ( \theta_\nu - \theta^0_\nu)$     (15)


$\displaystyle E^{str-bnd}_{\mu\nu} (r_\mu, \theta_\nu) = k^{sb}_{\mu\nu}
( r_\mu - r^0_\mu) ( \theta_\nu - \theta^0_\nu)$     (16)


$\displaystyle E^{str-Rtors}_{\mu\nu} (r_\mu, \phi_\nu) = k^{srt}_{\mu\nu}
( r_\mu - r^0_\mu) ( \phi_\nu - \phi^0_\nu)$     (17)

While similar expressions can be easily derived for all couplings involving stiff ICs, for the sake of clarity, we describe in more detail the effect of coupling a stiff coordinate with a soft one, as for instance a selected dihedral angle $ \delta^\prime$ and a neighboring stretching distance, $ r_1$. If we adopt a linear form for the term related to the stiff coordinate, $ f(r_1) = (r_1-r_1^0)$, and merge the coefficient and the function related to the soft coordinate into a generic function $ G(\delta\prime)$, equation (13) leads to

$\displaystyle E_{Coupl}\equiv E^{Coupl}(r_1,\delta^\prime) =
(r_1-r_1^{0})G(\delta^\prime)$     (18)

The shape of $ G(\delta\prime)$ should be flexible enough so as to reproduce the reference profiles obtained by QM calculations. In this sense, a Fourier series can again provide the required flexibility, together with other desirable features. [21] All requisites can be verified using the functional form:
$\displaystyle G(\delta^\prime) = \sum_i
k^{c}_i [1 + \sin(n_i \delta^\prime - \gamma^c_i)]$     (19)

which leads, to the following explicit coupling terms:
$\displaystyle E^{str-tors}_{\mu\nu}(r_\mu, \delta_\nu) = \sum^{N_{sin^i_\nu}}_{j=1}
k^{st}_{j\mu\nu} ( r_\mu - r^0_\mu)
\sin(n_j^\nu \delta_\nu - \gamma_j^\nu)$     (20)


$\displaystyle E^{bnd-tors}_{\mu\nu} (\theta_\mu, \delta_\nu) = \sum^{N_{sin^i_\...
..._{j\mu\nu} ( \theta_\mu - \theta^0_\mu)
\sin(n_j^\nu \delta_\nu - \gamma_j^\nu)$     (21)


$\displaystyle E^{Rtors-tors}_{\mu\nu} (\phi_\mu, \delta_\nu) = \sum^{N_{sin^i_\...
...rtt}_{j\mu\nu} ( \phi_\mu - \phi^0_\mu)
\sin(n_j^\nu \delta_\nu - \gamma_j^\nu)$     (22)


$\displaystyle E^{tors-tors}_{\mu, \nu}(\delta_\mu,\delta_\nu) =
\sum^{N_{sin^i_...
...tc} \sin(n_j^\mu \delta_\mu - \gamma_j^i)
\sin(m_k^\nu \delta_\nu - \gamma_k^i)$     (23)