generalizing potential energy penalties

pmetz

Hi All,

The built-in grs_interaction is handy for ionically bonded solids, but (assuming I understand it correctly) is not appropriate for molecular (covalently bonded) structures.

For e.g. an organic small molecule system with a variety of light atoms (C, N, S, H, ...), is there a simple way to control the coefficients for potentials defining different bonds and interactions? E.g. atomic_interaction will allow us to iterate over multiple sites, but I don't think there's a simple way to look up model coefficients depending on the identity of site1, site2, (...).

I can apply these constraints with foreknowledge of the system, i.e.

stretch_penalty(site1, site2, prm1, ...)

...

bend_penalty(site1, site2, site3, prm1, ...)

...

torsion_penalty(site1, site2, site3, site3, prm1, ...)

where prm1, ..., prmN could be tabulated e.g. in the universal force field potential set. This moves the modelling strategy away from the intersection of molecular dynamics and Rietveld methods, though, and requires lots of typing.

Thoughts or comments? Obvious things I've missed?

Thanks,
Peter

alancoelho

Hi Peter

Yes, grs_intereaction is not directional and hence not suited for covalent bonding.

I am not sure what your mean exactly however. Are you wanting to including torsion angles in penalties.

If yes then atomic_interaction can be used together with the AI_Cos_Angle and AI_Angle functions.

sites_geometry can also be used together with the functions:

Sites_Geometry_Distance($Name)
Sites_Geometry_Angle($Name)
Sites_Geometry_Dihedral_Angle($Name)

See examples:

test_examples\sites_geometry_1.inp
test_examples\sites_geometry_2.inp

If you get stuck then let me know, or, if you have ideas on how to improve these types of penalties then also let me know.

cheers
alan

pmetz

Hi Alan,

Thanks for the comments, and sorry for the late response.

I've abstracted the penalty creation part using sites_geometry and bond stretching and bending models based on [1] (macro examples below). The catch is that these function arguments are tabulated parameters for different elements, so as far as I know, the atomic interactions have to be defined a priori.

This is helpful for restraining molecule geometry, e.g., where we know well what connectivity is expected. The question I have is whether it's possible to look up the parameters for site1, site2, ... on the fly, so that bond penalties could be applied between structure fragments whose spatial relationships we don't know a priori.

E.g. hypothetical solution-- pseudocode using parameter lookup for bond stretches < cutoff radius:

for site2 in (sites < r_cutoff):
   stretch(site1,              <--- macro to insert bond stretch penalty
              site2,
              r_table(site1),  <--- equilibrium distance parameter
              r_table(site2),
              z_table(site1),  <--- effective charge parameter
              z_table(site2),
              ...
              )

For penalties defined a priori, I currently have the parameter values stored in a .inc file with a loader to populate the namespace of the input script:

macro load_prm(arg) {
    ''' invoke UFF prm instances for lab '''
    arg
}



' UFF parameters table from [Rappe 1992]
'
' ref: Rappe, A.K., Casewit, C.J., Colwell, K.S., Goddard III,
'        W.A., and Skiff, W.M. (1922) "UFF, a Full Periodic 
'        Table Force Field for Molecular Mechanics and 
'        Molecular Dynamics Simulations." J. Am. Chem. Soc.
'        114, 10024-10035.
'
' Table I. Atomic Data.
'
'              VALENCE                NONBOND          
'           ----------------  -----------------------  effect
'           bond    angle     dist    energy  scale    charge
' type      r_i     theta_i   x_i     D_i     Xi       Z_i
' -----------------------------------------------------------
macro uff_H      { prm !ri_H      0.3540     prm !thi_H      180.0000   prm !chii_H      0.89       prm !xi_H      2.8860     prm !di_H      0.0440     prm !si_H      12.0000    prm !zi_H      0.7120     }
...

Penalty example:
===========

macro stretch(si, sj, ri, rj, zi, zj, ci, cj, n, weight, dist0, dist, eval) {
    ''' bond stretch penalty (E_ij) '''
    local #m_unique !rij0 = get_rij( ri, rj, n, ci, cj ) ; : dist0
    sites_geometry #m_unique c1 load site_to_restrain { si sj }
        local #m_unique rij = Sites_Geometry_Distance(c1) ; : dist
        penalty = weight 
            IF rij > 2 rij0 THEN
                Er_harmonic(rij, zi, zj, ri, rj, n, ci, cj) ' large penalty at large distance
            ELSE
                Er_morse(rij, zi, zj, ri, rj, n, ci, cj) ' anharmonicity more accurate at short range
            ENDIF
            ; : eval
}

macro bend(si, sj, sk, th0, ri, rj, rk, zi, zk, ci, cj, ck, nij, njk, weight, ang, eval) {
    ''' bond bend penalty (E_ijk) '''
    sites_geometry #m_unique c1 load site_to_restrain { si sj sk }
        local #m_unique theta = Sites_Geometry_Angle(c1) ; : ang
        penalty = weight Etheta(theta, ri, rj, rk, th0, zi, zk, nij, njk, ci, cj, ck) ; : eval
}

Ref:
===
[1] A. K. Rappé, C. J. Casewit, K. S. Colwell, W. A. Goddard, and W. M. Skiff, “UFF, a Full Periodic Table Force Field for Molecular Mechanics and Molecular Dynamics Simulations,” J. Am. Chem. Soc., vol. 114, no. 25, pp. 10024–10035, 1992. https://pubs.acs.org/doi/abs/10.1021/ja00051a040