md.bond

Overview

md.bond.fene FENE bond potential.
md.bond.harmonic Harmonic bond potential.
md.bond.table Tabulated bond potential.

Details

Bond potentials.

Bonds add forces between specified pairs of particles and are typically used to model chemical bonds between two particles.

By themselves, bonds that have been specified in an initial configuration do nothing. Only when you specify an bond force (i.e. bond.harmonic), are forces actually calculated between the listed particles.

class hoomd.md.bond.coeff

Define bond coefficients.

The coefficients for all bond potentials are specified using this class. Coefficients are specified per bond type.

There are two ways to set the coefficients for a particular bond potential. The first way is to save the bond potential in a variable and call set() directly. See below for an example of this.

The second method is to build the coeff class first and then assign it to the bond potential. There are some advantages to this method in that you could specify a complicated set of bond potential coefficients in a separate python file and import it into your job script.

Example:

my_coeffs = hoomd.md.bond.coeff();
my_bond_force.bond_coeff.set('polymer', k=330.0, r=0.84)
my_bond_force.bond_coeff.set('backbone', k=330.0, r=0.84)
set(type, **coeffs)

Sets parameters for bond types.

Parameters:
  • type (str) – Type of bond (or a list of type names)
  • coeffs – Named coefficients (see below for examples)

Calling set() results in one or more parameters being set for a bond type. Types are identified by name, and parameters are also added by name. Which parameters you need to specify depends on the bond potential you are setting these coefficients for, see the corresponding documentation.

All possible bond types as defined in the simulation box must be specified before executing run(). You will receive an error if you fail to do so. It is not an error, however, to specify coefficients for bond types that do not exist in the simulation. This can be useful in defining a potential field for many different types of bonds even when some simulations only include a subset.

Examples:

my_bond_force.bond_coeff.set('polymer', k=330.0, r0=0.84)
my_bond_force.bond_coeff.set('backbone', k=1000.0, r0=1.0)
my_bond_force.bond_coeff.set(['bondA','bondB'], k=100, r0=0.0)

Note

Single parameters can be updated. If both k and r0 have already been set for a particle type, then executing coeff.set('polymer', r0=1.0) will update the value of r0 and leave the other parameters as they were previously set.

class hoomd.md.bond.fene(name=None)

FENE bond potential.

Parameters:name (str) – Name of the bond instance.

fene specifies a FENE potential energy between the two particles in each defined bond.

\[V(r) = - \frac{1}{2} k r_0^2 \ln \left( 1 - \left( \frac{r - \Delta}{r_0} \right)^2 \right) + V_{\mathrm{WCA}}(r)\]

where \(\vec{r}\) is the vector pointing from one particle to the other in the bond, \(\Delta = (d_i + d_j)/2 - 1\), \(d_i\) is the diameter of particle \(i\), and

\begin{eqnarray*} V_{\mathrm{WCA}}(r) = & 4 \varepsilon \left[ \left( \frac{\sigma}{r - \Delta} \right)^{12} - \left( \frac{\sigma}{r - \Delta} \right)^{6} \right] + \varepsilon & r-\Delta < 2^{\frac{1}{6}}\sigma\\ = & 0 & r-\Delta \ge 2^{\frac{1}{6}}\sigma \end{eqnarray*}

Coefficients:

  • \(k\) - attractive force strength k (in units of energy/distance^2)
  • \(r_0\) - size parameter r0 (in distance units)
  • \(\varepsilon\) - repulsive force strength epsilon (in energy units)
  • \(\sigma\) - repulsive force interaction distance sigma (in distance units)

Examples:

fene = bond.fene()
fene.bond_coeff.set('polymer', k=30.0, r0=1.5, sigma=1.0, epsilon= 2.0)
fene.bond_coeff.set('backbone', k=100.0, r0=1.0, sigma=1.0, epsilon= 2.0)
disable(log=False)

Disable the force.

Parameters:log (bool) – Set to True if you plan to continue logging the potential energy associated with this force.

Examples:

force.disable()
force.disable(log=True)

Executing the disable command will remove the force from the simulation. Any hoomd.run() command executed after disabling a force will not calculate or use the force during the simulation. A disabled force can be re-enabled with enable().

By setting log to True, the values of the force can be logged even though the forces are not applied in the simulation. For forces that use cutoff radii, setting log=True will cause the correct r_cut values to be used throughout the simulation, and therefore possibly drive the neighbor list size larger than it otherwise would be. If log is left False, the potential energy associated with this force will not be available for logging.

enable()

Enable the force.

Examples:

force.enable()

See disable().

get_energy(group)

Get the energy of a particle group.

Parameters:group (hoomd.group) – The particle group to query the energy for.
Returns:The last computed energy for the members in the group.

Examples:

g = group.all()
energy = force.get_energy(g)
get_net_force(group)

Get the force of a particle group.

Parameters:group (hoomd.group) – The particle group to query the force for.
Returns:The last computed force for the members in the group.

Examples

g = group.all() force = force.get_net_force(g)

get_net_virial(group)

Get the virial of a particle group.

Parameters:group (hoomd.group) – The particle group to query the virial for.
Returns:The last computed virial for the members in the group.

Examples

g = group.all() virial = force.get_net_virial(g)

class hoomd.md.bond.harmonic(name=None)

Harmonic bond potential.

Parameters:name (str) – Name of the bond instance.

harmonic specifies a harmonic potential energy between the two particles in each defined bond.

\[V(r) = \frac{1}{2} k \left( r - r_0 \right)^2\]

where \(\vec{r}\) is the vector pointing from one particle to the other in the bond.

Coefficients:

  • \(k\) - force constant k (in units of energy/distance^2)
  • \(r_0\) - bond rest length r0 (in distance units)

Example:

harmonic = bond.harmonic(name="mybond")
harmonic.bond_coeff.set('polymer', k=330.0, r0=0.84)
disable(log=False)

Disable the force.

Parameters:log (bool) – Set to True if you plan to continue logging the potential energy associated with this force.

Examples:

force.disable()
force.disable(log=True)

Executing the disable command will remove the force from the simulation. Any hoomd.run() command executed after disabling a force will not calculate or use the force during the simulation. A disabled force can be re-enabled with enable().

By setting log to True, the values of the force can be logged even though the forces are not applied in the simulation. For forces that use cutoff radii, setting log=True will cause the correct r_cut values to be used throughout the simulation, and therefore possibly drive the neighbor list size larger than it otherwise would be. If log is left False, the potential energy associated with this force will not be available for logging.

enable()

Enable the force.

Examples:

force.enable()

See disable().

get_energy(group)

Get the energy of a particle group.

Parameters:group (hoomd.group) – The particle group to query the energy for.
Returns:The last computed energy for the members in the group.

Examples:

g = group.all()
energy = force.get_energy(g)
get_net_force(group)

Get the force of a particle group.

Parameters:group (hoomd.group) – The particle group to query the force for.
Returns:The last computed force for the members in the group.

Examples

g = group.all() force = force.get_net_force(g)

get_net_virial(group)

Get the virial of a particle group.

Parameters:group (hoomd.group) – The particle group to query the virial for.
Returns:The last computed virial for the members in the group.

Examples

g = group.all() virial = force.get_net_virial(g)

class hoomd.md.bond.table(width, name=None)

Tabulated bond potential.

Parameters:
  • width (int) – Number of points to use to interpolate V and F
  • name (str) – Name of the potential instance

table specifies that a tabulated bond potential should be applied between the two particles in each defined bond.

The force \(\vec{F}\) is (in force units) and the potential \(V(r)\) is (in energy units):

\begin{eqnarray*} \vec{F}(\vec{r}) = & 0 & r < r_{\mathrm{min}} \\ = & F_{\mathrm{user}}(r)\hat{r} & r < r_{\mathrm{max}} \\ = & 0 & r \ge r_{\mathrm{max}} \\ \\ V(r) = & 0 & r < r_{\mathrm{min}} \\ = & V_{\mathrm{user}}(r) & r < r_{\mathrm{max}} \\ = & 0 & r \ge r_{\mathrm{max}} \\ \end{eqnarray*}

where \(\vec{r}\) is the vector pointing from one particle to the other in the bond. Care should be taken to define the range of the bond so that it is not possible for the distance between two bonded particles to be outside the specified range. On the CPU, this will throw an error. On the GPU, this will throw an error if GPU error checking is enabled.

\(F_{\mathrm{user}}(r)\) and \(V_{\mathrm{user}}(r)\) are evaluated on width grid points between \(r_{\mathrm{min}}\) and \(r_{\mathrm{max}}\). Values are interpolated linearly between grid points. For correctness, you must specify the force defined by: \(F = -\frac{\partial V}{\partial r}\)

The following coefficients must be set for each bond type:

  • \(F_{\mathrm{user}}(r)\) and \(V_{\mathrm{user}}(r)\) - evaluated by func (see example)
  • coefficients passed to func - coeff (see example)
  • \(r_{\mathrm{min}}\) - rmin (in distance units)
  • \(r_{\mathrm{max}}\) - rmax (in distance units)

The table width is set once when bond.table is specified. There are two ways to specify the other parameters.

Set table from a given function

When you have a functional form for V and F, you can enter that directly into python. table will evaluate the given function over width points between rmin and rmax and use the resulting values in the table:

def harmonic(r, rmin, rmax, kappa, r0):
   V = 0.5 * kappa * (r-r0)**2;
   F = -kappa*(r-r0);
   return (V, F)

btable = bond.table(width=1000)
btable.bond_coeff.set('bond1', func=harmonic, rmin=0.2, rmax=5.0, coeff=dict(kappa=330, r0=0.84))
btable.bond_coeff.set('bond2', func=harmonic, rmin=0.2, rmax=5.0, coeff=dict(kappa=30, r0=1.0))

Set a table from a file

When you have no function for for V or F, or you otherwise have the data listed in a file, table can use the given values directly. You must first specify the number of rows in your tables when initializing bond.table. Then use set_from_file() to read the file:

btable = bond.table(width=1000)
btable.set_from_file('polymer', 'btable.file')

Note

For potentials that diverge near r=0, make sure to set rmin to a reasonable value. If a potential does not diverge near r=0, then a setting of rmin=0 is valid.

Note

Ensure that rmin and rmax cover the range of possible bond lengths. When gpu error checking is on, a error will be thrown if a bond distance is outside than this range.

disable(log=False)

Disable the force.

Parameters:log (bool) – Set to True if you plan to continue logging the potential energy associated with this force.

Examples:

force.disable()
force.disable(log=True)

Executing the disable command will remove the force from the simulation. Any hoomd.run() command executed after disabling a force will not calculate or use the force during the simulation. A disabled force can be re-enabled with enable().

By setting log to True, the values of the force can be logged even though the forces are not applied in the simulation. For forces that use cutoff radii, setting log=True will cause the correct r_cut values to be used throughout the simulation, and therefore possibly drive the neighbor list size larger than it otherwise would be. If log is left False, the potential energy associated with this force will not be available for logging.

enable()

Enable the force.

Examples:

force.enable()

See disable().

get_energy(group)

Get the energy of a particle group.

Parameters:group (hoomd.group) – The particle group to query the energy for.
Returns:The last computed energy for the members in the group.

Examples:

g = group.all()
energy = force.get_energy(g)
get_net_force(group)

Get the force of a particle group.

Parameters:group (hoomd.group) – The particle group to query the force for.
Returns:The last computed force for the members in the group.

Examples

g = group.all() force = force.get_net_force(g)

get_net_virial(group)

Get the virial of a particle group.

Parameters:group (hoomd.group) – The particle group to query the virial for.
Returns:The last computed virial for the members in the group.

Examples

g = group.all() virial = force.get_net_virial(g)

set_from_file(bondname, filename)

Set a bond pair interaction from a file.

Parameters:
  • bondname (str) – Name of bond
  • filename (str) – Name of the file to read

The provided file specifies V and F at equally spaced r values. Example:

#r  V    F
1.0 2.0 -3.0
1.1 3.0 -4.0
1.2 2.0 -3.0
1.3 1.0 -2.0
1.4 0.0 -1.0
1.5 -1.0 0.0

The first r value sets rmin, the last sets rmax. Any line with # as the first non-whitespace character is is treated as a comment. The r values must monotonically increase and be equally spaced. The table is read directly into the grid points used to evaluate \(F_{\mathrm{user}}(r)\) and \(V_{\mathrm{user}}(r)\).