Non-bonded interactions¶

Overview

The pygamd MD engine provides a few of functions for non-bonded interactions. However, it supports well self-defined analytical functions via writting codes of device function in script.

Non-bonded functions	`force.nonbonded`
Self-defined functions	`force.nonbonded`
Non-bonded functions for charged beads	`force.nonbonded_c`
Self-defined functions for charged beads	`force.nonbonded_c`

Non-bonded functions¶

Description:

The function describing non-bonded interactions could be either the one called from non-bonded interaction function libary, or the one defined by user himself in script. Non-bonded interaction function libary contains Lennard-Jones function named as ‘lj’ and harmonic function named as ‘harmonic’.

Lennard-Jones function (‘lj’)
\begin{eqnarray*} V_{\mathrm{LJ}}(r) = & 4 \epsilon \left[ \left( \frac{\sigma}{r} \right)^{12} - \alpha \left( \frac{\sigma}{r} \right)^{6} \right] & r < r_{\mathrm{cut}} \\ = & 0 & r \ge r_{\mathrm{cut}} \\ \end{eqnarray*}
The following coefficients must be set per unique pair of particle types:

\(\epsilon\) - the depth of the potential well (in energy units)

\(\sigma\) - the collision diameter (in distance units)

\(\alpha\) - the factor of attraction (unitless)

\(r_{\mathrm{cut}}\) - cutoff radius (in distance units) - note: equal to or smaller than the global rcut specified in force.nonbonded

Harmonic function (‘harmonic’)
\begin{eqnarray*} V_{\mathrm{H}}(r) = \frac{1}{2}\alpha \left( r - r_{\mathrm{cut}} \right)^{2}, r < r_{\mathrm{cut}} \end{eqnarray*}
Coefficients:

\(\alpha\) - spring constant (in units of energy/distance^2)

\(r_{\mathrm{cut}}\) - cutoff radius (in distance units) - note: equal to or smaller than the global rcut specified in force.nonbonded

class force.nonbonded(info, rcut, func, exclusion=None)¶

Constructor of non-bonded interaction calculation object.

Parameters:

info – system information.
rcut – cut-off radius of interactions.
func – function name.
exclusion – a python list of exclusions, the candidates are ‘bond’, ‘angle’, ‘dihedral’, the default is None.

setParams(type_i, type_j, param)¶: specifies interaction parameters with type_i, type_j, a list of parameters.

Example:

fn = pygamd.force.nonbonded(info=mst, rcut=3.0, func='lj')
fn.setParams(type_i="a", type_j="a", param=[1.0, 1.0, 1.0, 3.0])
app.add(fn)


fn = pygamd.force.nonbonded(info=mst, rcut=3.0, func='lj', exclusion=['bond'])
fn.setParams(type_i="a", type_j="a", param=[1.0, 1.0, 1.0, 3.0])
app.add(fn)

Self-defined functions¶

Description:

The device function for non-bonded interactions could be written in script and conveyed to kernel funciton for calculation. The function has three parameters where rsq, param, and fp are square of distance, interaction parameters, and force and potential, respectively.

With the potential form of non-bonded interactions \(p(r)\), the expression of parameters in script are:

p = \(p(r)\)
f = \(-(\partial p(r)/\partial r)(1/r)\)

Function code template:

@cuda.jit(device=True)
def func(rsq, param, fp):
        rcut = param[0]
        p1 = param[1]
        p2 = param[2]
        p3 = param[3]
        ...
        if rsq<rcut*rcut:
                calculation codes
                ...
                fp[0]=f
                fp[1]=p

fn = pygamd.force.nonbonded(info, rcut, func)
fn.setParams(type_i, type_j, param=[rcut, p1, p2, p3, ...])
....
app.add(fn)

Example:

from numba import cuda
import numba as nb

@cuda.jit(device=True)
def lj(rsq, param, fp):
        epsilon = param[0]
        sigma = param[1]
        alpha = param[2]
        rcut = param[3]
        if rsq<rcut*rcut:
                sigma2 = sigma*sigma
                r2inv = sigma2/rsq;
                r6inv = r2inv * r2inv * r2inv;
                f = nb.float32(4.0) * epsilon * r2inv * r6inv * (nb.float32(12.0)
                        * r6inv - nb.float32(6.0) * alpha)/sigma2
                p = nb.float32(4.0) * epsilon * r6inv * ( r6inv - nb.float32(1.0))
                fp[0]=f
                fp[1]=p

fn = pygamd.force.nonbonded(info=mst, rcut=3.0, func=lj)
fn.setParams(type_i="a", type_j="a", param=[1.0, 1.0, 1.0, 3.0])
app.add(fn)

Non-bonded functions for charged beads¶

Description:

The function describing non-bonded interactions including electrostatic part could be either the one called from non-bonded interaction function libary, or the one defined by user himself in script.

Lennard-Jones-coulomb function (‘lj_coulomb’)
\begin{eqnarray*} V_{\mathrm{LJ}}(r) = & 4 \epsilon \left[ \left( \frac{\sigma}{r} \right)^{12} - \alpha \left( \frac{\sigma}{r} \right)^{6} \right] + f\frac{q_{i} q_{j}}{\epsilon_{r}r} & r < r_{\mathrm{cut}} \\ = & 0 & r \ge r_{\mathrm{cut}} \\ \end{eqnarray*}
Following coefficients need being set per pair of particle types:

\(\epsilon\) - the depth of the potential well (in energy units)

\(\sigma\) - the collision diameter (in distance units)

\(\alpha\) - the factor of attraction (unitless)

\(\epsilon_{r}\) - dielectric coefficient : \(f= 1/4\pi \epsilon_0=138.935\text{ }kJ\text{ }mol^{-1}\text{ }nm\text{ }e^{-2}\)

\(r_{\mathrm{cut}}\) - cutoff radius (in distance units) - note: equal to or smaller than the global rcut specified in force.nonbonded_c

class force.nonbonded_c(info, rcut, func, exclusion=None)¶

Constructor of non-bonded interaction calculation object.

Parameters:

info – system information.
rcut – cut-off radius of interactions.
func – function name.
exclusion – a python list of exclusions, the candidates are ‘bond’, ‘angle’, ‘dihedral’, the default is None.

setParams(type_i, type_j, param)¶: specifies interaction parameters with type_i, type_j, a list of parameters.

Example:

fn = pygamd.force.nonbonded_c(info=mst, rcut=3.0, func='lj_coulomb')
fn.setParams(type_i="a", type_j="a", param=[1.0, 1.0, 1.0, 15.0, 3.0])
app.add(fn)


fn = pygamd.force.nonbonded_c(info=mst, rcut=3.0, func='lj_coulomb', exclusion=['bond'])
fn.setParams(type_i="a", type_j="a", param=[1.0, 1.0, 1.0, 15.0, 3.0])
app.add(fn)

Self-defined functions for charged beads¶

Description:

The device function for non-bonded interactions could be written in script and conveyed to kernel funciton for calculation. The function has three parameters where rsq, qi, qj, param, and fp are the square of distance, qi and qj, interaction parameters, and force and potential, respectively.

With the potential form of non-bonded interactions \(p(r)\), the expression of parameters in script are:

p = \(p(r)\)
f = \(-(\triangle p(r)/\triangle r)(1/r)\)

Function code template:

@cuda.jit(device=True)
def func(rsq, qi, qj, param, fp):
        rcut = param[0]
        p1 = param[1]
        p2 = param[2]
        p3 = param[3]
        ...
        if rsq<rcut*rcut:
                calculation codes
                ...
                fp[0]=f
                fp[1]=p

fn = pygamd.force.nonbonded_c(info, rcut, func)
fn.setParams(type_i, type_j, param=[rcut, p1, p2, p3, ...])
....
app.add(fn)

Example:

from numba import cuda
import numba as nb

@cuda.jit(device=True)
def lj_coulomb(rsq, qi, qj, param, fp):
        epsilon = param[0]
        sigma = param[1]
        alpha = param[2]
        epsilonr = param[3]
        rcut = param[4]
        coulomb_eff = 138.935/epsilonr
        if rsq<rcut*rcut:
                sigma2 = sigma*sigma
                r2inv = sigma2/rsq;
                r6inv = r2inv * r2inv * r2inv;
                f = nb.float32(4.0) * epsilon * r2inv * r6inv * (nb.float32(12.0)
                        * r6inv - nb.float32(6.0) * alpha)/sigma2 + coulomb_eff*qi*qj*r2inv*rinv
                p = nb.float32(4.0) * epsilon * r6inv * ( r6inv - nb.float32(1.0))
                        + coulomb_eff*qi*qj*rinv
                fp[0]=f
                fp[1]=p

fn = pygamd.force.nonbonded_c(info=mst, rcut=3.0, func=lj_coulomb)
fn.setParams(type_i="a", type_j="a", param=[1.0, 1.0, 1.0, 15.0, 3.0])
app.add(fn)