DreidingPotentialBuilder¶
- class DreidingPotentialBuilder(use_bonds=None, include_bond_stretching=None, include_bond_angles=None, include_torsions=None, include_inversions=None, include_h_bonds=None, include_lennard_jones=None, include_electrostatic=None, lennard_jones_cutoff=None, lennard_jones_smoothing_length=None, h_bond_scale=None, h_bond_cutoff=None, electrostatic_pairwise_cutoff=None, electrostatic_smoothing_length=None, electrostatic_spme_cutoff=None, electrostatic_accuracy=None, electrostatic_scale=None, use_electrostatic_direct_sum=None, ignore_missing_types=None)¶
This object builds a DREIDING potential using TremoloX.
- Parameters:
use_bonds – Controls if any bonding (stretching, angle, torsion and inversion) terms are to be included. Default: True.
include_bond_stretching (bool) – Controls if bonding stretching terms are included. Default: use_bonds.
include_bond_angles (bool) – Controls if angle bending terms are included. Default: use_bonds.
include_torsions (bool) – Controls if bond torsion terms are included. Default: use_bonds.
include_inversions (bool) – Controls if bond inversion torsion terms are included. Default: use_bonds.
include_h_bonds (bool) – Controls if hydrogen bonding terms are included. Default: use_bonds.
include_lennard_jones (bool) – Controls if Lennard-Jones terms are included. Default: True.
include_electrostatic (bool) – Controls if electrostatic terms are included. Default: True.
lennard_jones_cutoff (PhysicalQuantity of type distance) – Distance beyond which the Lennard-Jones terms are truncated. Default: 10.0 * Angstrom.
lennard_jones_smoothing_length (PhysicalQuantity of type distance) – Distance over which the Lennard-Jones terms are brought to zero at the cutoff. Default: 2.0 * Angstrom.
h_bond_scale (float) – Scaling for the strength of the hydrogen bond term to allow for setting atomic charges. Default: 1.0.
h_bond_cutoff (PhysicalQuantity of type distance) – Distance beyond which the Hydrogen bond terms are truncated. Default: 7.5 * Angstrom.
electrostatic_pairwise_cutoff (PhysicalQuantity of type distance) – Distance at which the pairwise electrostatic terms are brought to zero in pairwise summation. Default: 12 * Angstrom
electrostatic_smoothing_length (PhysicalQuantity of type distance) – Distance over which the electrostatic terms are brought to zero at the cutoff when using the direct electrostatic sum. Default: 2.0 * Angstrom
electrostatic_spme_cutoff (PhysicalQuantity of type distance) – Distance at which the pairwise electrostatic terms are brought to zero in SPME summation. Default: 7.5 * Angstrom
electrostatic_accuracy (float) – Sets the relative accuracy of the SPME summation. Default: 0.0001
electrostatic_scale (float) – Set the relative strength of the electrostatic terms. Default: 1.0
use_electrostatic_direct_sum (bool) – Use the direct Coulomb sum to calculate the electrostatic terms in MoleculeConfigurations. Default: False
ignore_missing_types (bool) – Whether or not molecules missing Dreiding types are ignored. Default: False
- assignAtomTypes(configuration, overwrite=True)¶
Function that gets the atom types based on the bonding pattern.
- Parameters:
configuration (AtomicConfiguration) – The configuration of the system that is having types assigned. Default: True
overwrite (bool) – Whether or not the existing Dreiding types are overwritten.
- checkAtomTypes(configuration, add_missing_terms=False)¶
Function that tests that each atom has only one valid Dreiding type.
- Parameters:
configuration (AtomicConfiguration) – The configuration whose types are being checked.
add_missing_terms (bool) – Whether or not to add terms for missing Dreiding types. Default: False
- createCalculator(configuration, check_types=True, atomic_charges=None, charge_averaging=None)¶
Create and return a calculator with the Dreiding potential.
- Parameters:
configuration (AtomicConfiguration) – The configuration the potential is to be used with.
check_types (bool) – Flag as to whether or not the correctness of the types are checked. Default: True
atomic_charges (PhysicalQuantity of type charge) – Custom atomic charges to be used in the potential.
charge_averaging (NLFlag | None) – The type of averaging of charges to be used, if any. Accepted values are None, AtomConnectivity to average over atoms bonded to the same elements and AtomTypes to average over atoms of the same forcefield type. Default: None
- Returns:
The calculator containing the potential.
- Return type:
- dreidingType(index, configuration, add_missing_terms=False)¶
Gives the Dreiding type on an atom in a configuration, ignoring other tags.
- Parameters:
index (int) – Index of the atom for which the type is sought.
configuration (AtomicConfiguration) – The configuration containing the atom.
add_missing_terms (bool) – Whether or not generic Dreiding types are considered.
- Returns:
The type of the atom.
- Return type:
str
- electrostaticAccuracy()¶
- Returns:
The relative accuracy of the electrostatic summation for SPME.
- Return type:
float
- electrostaticPairwiseCutoff()¶
- Returns:
The cutoff used in pairwise electrostatic summation.
- Return type:
PhysicalQuantity (length)
- electrostaticScale()¶
- Returns:
Scaling factor for the electrostatic terms.
- Return type:
float
- electrostaticSmoothingLength()¶
- Returns:
Distance at which electrostatic terms are splined when using direct summation.
- Return type:
PhysicalQuantity (length)
- electrostaticSpmeCutoff()¶
- Returns:
The cutoff used in SPME electrostatic summation.
- Return type:
PhysicalQuantity (length)
- hBondCutoff()¶
- Returns:
Distance at which hydrogen bond terms are ignored.
- Return type:
PhysicalQuantity (length)
- hBondScale()¶
- Returns:
Scaling of the hydrogen bond energy because of included charges.
- Return type:
float
- ignoreMissingTypes()¶
- Returns:
Whether or not missing types are ignored on molecules.
- Return type:
bool
- includeBondAngles()¶
- Returns:
Whether or not bond angle terms are included in the potential set.
- Return type:
bool
- includeBondStretching()¶
- Returns:
Whether or not bond stretching terms are included in the potential set.
- Return type:
bool
- includeElectrostatic()¶
- Returns:
Whether or not electrostatic terms are included in the potential set.
- Return type:
bool
- includeHBonds()¶
- Returns:
Whether or not hydrogen-bond terms are included in the potential set.
- Return type:
bool
- includeInversions()¶
- Returns:
Whether or not bond inversion terms are included in the potential set.
- Return type:
bool
- includeLennardJones()¶
- Returns:
Whether or not Lennard Jones terms are included in the potential set.
- Return type:
bool
- includeTorsions()¶
- Returns:
Whether or not bond torsion terms are included in the potential set.
- Return type:
bool
- isDreidingType(type_label)¶
Function that tests if a type is in the potential
- Parameters:
type_label (str) – The type being tested.
- Returns:
Whether or not the tag is defined in the potential.
- Return type:
bool
- lennardJonesCutoff()¶
- Returns:
Distance at which Lennard Jones terms are ignored.
- Return type:
PhysicalQuantity (length)
- lennardJonesSmoothingLength()¶
- Returns:
Distance at which Lennard Jones terms are splined.
- Return type:
PhysicalQuantity (length)
- uniqueString()¶
Return a unique string representing the state of the object.
- useBonds()¶
- Returns:
Whether or not bond terms are included in the potential set.
- Return type:
bool
- useElectrostaticDirectSum()¶
- Returns:
Whether or not direct Coulomb sums are used for molecules.
- Return type:
bool
Usage Examples¶
Investigate the diffusion of chloroform molecules in LTA zeolites using the DREIDING potential.
# -------------------------------------------------------------
# Bulk Configuration
# -------------------------------------------------------------
bulk_configuration = nlread('LTA_Chloroform.hdf5')[-1]
# -------------------------------------------------------------
# DREIDING Calculator
# -------------------------------------------------------------
potential_builder = DreidingPotentialBuilder()
potential_builder.assignAtomTypes(bulk_configuration)
calculator = potential_builder.createCalculator(bulk_configuration)
bulk_configuration.setCalculator(calculator)
# -------------------------------------------------------------
# Molecular Dynamics
# -------------------------------------------------------------
initial_velocity = MaxwellBoltzmannDistribution(
temperature=300.0*Kelvin,
remove_center_of_mass_momentum=True
)
method = NVTNoseHoover(
time_step=1*femtoSecond,
reservoir_temperature=300*Kelvin,
thermostat_timescale=100*femtoSecond,
heating_rate=0*Kelvin/picoSecond,
chain_length=3,
initial_velocity=initial_velocity,
)
constraints = [FixCenterOfMass()]
md_trajectory = MolecularDynamics(
bulk_configuration,
constraints=constraints,
trajectory_filename='LTA_Chloroform_Trajectory.hdf5',
steps=100000,
log_interval=500,
method=method
)
bulk_configuration = md_trajectory.lastImage()
Notes¶
The DreidingPotentialBuilder class enables the creation of a TremoloXCalculator
that
implements the DREIDING potential of Mayo, Olafson and Goddard[1]. This is a bonded
valence forcefield that represents the potential energy of a molecule as a collection of simple
functions based on bond lengths, bond angles, torsion angles, inversion angles and inter-atomic
distances. The overall mathematical form of the DREIDING potential can be given as:
Here \(r\), \(\theta\), \(\phi\) and \(\chi\) represent inter-atomic and bond distances, angles, torsions and inversions respectively in the configuration.
To create a DREIDING calculator, an instance of the DreidingPotentialBuilder
must first
be created. The constructor for this class has a number of options that controls the way in which
the calculator is built. Individual types of interactions can be excluded by setting the include_
options to False
. By default all possible terms are included. All valence terms can also be
excluded by setting the argument use_bonds
to False. This is appropriate for configurations
consisting of single atoms that are not bonded to each other.
As the DreidingPotentialBuilder
class uses the bonding in the configuration to determine
appropriate potentials, bonds should be set on the configuration. This can be done with the
findBonds()
method in either the MoleculeConfiguration
or
BulkConfiguration
class. This automatically assigns bonds between atoms that are
within the covalent radii of each other. Specific bonds can also be assigned with the
setBonds()
, which takes a list of the indices of pairs of atoms that are
to be considered bonded.
There are then a number of options that control the non-bonding interactions. For hydrogen bonding,
the cutoff used in calculating hydrogen bonding interactions can be set with the parameter
h_bond_cutoff
. The strength of the hydrogen bonding interactions can also be scaled using the
h_bond_scale
parameter. In the original DREIDING paper it was mentioned that it may be
beneficial to scale down the strength of hydrogen bonding when combined with atomic partial charges,
as part of the hydrogen bonding effect may be replicated by the electrostatic interaction. By
default, the DREIDING potential builder adds unscaled hydrogen bonding potentials to the calculator.
Lennard-Jones interactions can be similarly controlled. The parameter lennard_jones_cutoff
sets
the distance beyond which Lennard-Jones interactions are truncated. To avoid integration problems
in molecular dynamics, Lennard-Jones interactions are brought to zero at the cutoff using a spline
function. The width of this spline can be set using the lennard_jones_smoothing_length
parameter.
Electrostatic interactions are a little more complicated. By default the
DreidingPotentialBuilder
adds a CoulombSPME
electrostatic solver for
calculators to be used with a bulk configuration, and a CoulombDSF
electrostatic solver
to calculators to be used with a molecule configuration. When adding an CoulombSPME
solver the cutoff used for calculating the real-space interactions can be set using the
electrostatic_spme_cutoff
argument. Likewise the accuracy of the smoothed particle mesh Ewald
solver can be set using the electrostatic_accuracy
argument. In the case where a
CoulombDSF
solver is added to the potential, the cutoff used can be set using the
electrostatic_pairwise_cutoff
argument. In cases where molecular electrostatics are to be
compared with values calculated in a bulk configuration, an explicit splined Coulomb potential
can be used for molecules instead of the damped shifted force potential. In this case a
CoulombN2Spline
electrostatic solver can be added by using
use_electrostatic_direct_sum
argument. As with the damped shifted force potential the cutoff is
controlled using the electrostatic_pairwise_cutoff
argument, and like the Lennard-Jones
potential is brought to zero with a spline whose length is defined by the
electrostatic_smoothing_length
argument. For all Coulomb solvers the total electrostatic
energies and forces can be scaled up or down using the electrostatic_scale
argument.
Once the DreidingPotentialBuilder
object is instantiated, a calculator is created using
the createCalculator
function. This takes a configuration and returns the appropriate DREIDING
calculator for that configuration. To use electrostatic potentials in the calculator, atomic
partial charges have to be assigned to each atom. As the DREIDING forcefield does not define
specific charges, by default charges are calculated with the QEqAtomicCharges
class.
This implements the QEq charge equilibration method of Rappe and Goddard[2].
These charges are then added to the calculator. It is also possible to add pre-calculated charges
from different sources using the atomic_charges
argument. This argument takes a list of charges
that have the same ordering as the atoms in the configuration. The partial atomic charges can also
be averaged according to a few different schemes. This is controlled by the charge_averaging
parameter, which takes an NLFlag which indicates the type of averaging. The flag
AtomConnectivity
averages charges on each atom based on the element and the other elements it is
bonded to. The flag AtomType
averages according to the assigned DREIDING type.
When creating a DREIDING calculator for a configuration, the configuration must have the
appropriate tags for each atom type. Tags used by the DREIDING forcefield all begin with the prefix
DREI_
. The next two characters in the tag are for the atomic symbol, followed by a possible
character to indicate the hybridization of the atom. Tags can be assigned manually, or they can be
assigned automatically. Automatic typing can be done with the assignAtomTypes
function. This
takes the configuration to have types added, as well as a possible flag to overwrite existing types.
The atom types are determined by the element and in some cases the bonding of the atom. The
validity of the existing types on a configuration can be checked by calling the checkAtomTypes
function. This function raises a NLValueError
if a problem is found with the tags on the
given configuration. The DREIDING type of a specific atom can also be found using the
dreidingType
function.