ReaxFFAtomicCharges

class ReaxFFAtomicCharges(r_cut=PhysicalQuantity(10.0, Ang), chi=None, eta=None, shielding=None)

This object calculates atomic charges based on the ReaxFF charge equilibration method.

Parameters:
  • r_cut (PhysicalQuantity of type length) – The cutoff used in determining charge-charge interactions.
    Default: 10 Angstrom.

  • chi (dict) – Custom atomic electronegativies, as a dictionary with PeriodicTableElement objects as keys and values as PhysicalQuantity of type energy.
    Default: Electronegativities from the ReaxFF CHONSSiCaCsKSrNaMgAlCu_2015 potential set.

  • eta (dict) – Custom atomic hardnesses, as a dictionary with PeriodicTableElement objects as keys and values as PhysicalQuantity of type energy.
    Default: Atomic hardnesses from the ReaxFF CHONSSiCaCsKSrNaMgAlCu_2015 potential set.

  • shielding (dict) – Custom electrostatic shielding parameters, as a dictionary with PeriodicTableElement objects as keys and values as PhysicalQuantity of type inverse distance.
    Default: Shielding parameters from the ReaxFF CHONSSiCaCsKSrNaMgAlCu_2015 potential set.

charges(configuration, total_charge=None, decompose_bulk=True)

Calculate the atomic charges on a given configuration.

Parameters:
  • configuration (AtomicConfiguration) – The configuration whose atomic charges are being calculated.

  • total_charge (PhysicalQuantity of type charge) – The total charge of the configuration.
    Default: Match the integer charge of the configuration.

  • decompose_bulk (bool) – Whether or not the configuration, if a BulkConfiguration, is decomposed into its parts to obtain the charges.
    Default: True

Returns:

An array of atomic charges.

Return type:

PhysicalQuantity of type charge

Usage Examples

Define a simulation box of 900 methanol molecules with a Dreiding calculator using ReaxFF calculated charges


# --------------------------------------------------------------------------------------------------
# Create Methanol Molecule
# --------------------------------------------------------------------------------------------------
# Define elements
elements = [Carbon, Hydrogen, Hydrogen, Hydrogen, Oxygen, Hydrogen]

# Define coordinates
cartesian_coordinates = [[ 0.048333359398, -0.11839206808 , -0.068353692376],
                         [ 1.138333281204, -0.11839206808 , -0.068353692376],
                         [-0.31499994787 , -0.11839206808 , -1.096015473978],
                         [-0.31499994787 , -1.008373277446,  0.445477198425],
                         [-0.363333307269,  0.889981209366,  0.513830890801],
                         [-0.013333332377,  1.747302557837,  0.018856179479]]*Angstrom

# Set up configuration
methanol = MoleculeConfiguration(
    elements=elements,
    cartesian_coordinates=cartesian_coordinates
)

# Add tags
methanol.addTags('DREI_C_3',  [0])
methanol.addTags('DREI_H_',   [1, 2, 3])
methanol.addTags('DREI_H__A', [5])
methanol.addTags('DREI_O_3',  [4])

# Add bonds
bonds = [[0, 1],
         [0, 2],
         [0, 3],
         [0, 4],
         [4, 5]]
methanol.setBonds(bonds)

# --------------------------------------------------------------------------------------------------
# Calculate Atomic Charges
# --------------------------------------------------------------------------------------------------
# Create the charger and calculate the charges
charger = ReaxFFAtomicCharges()
methanol_charges = charger.charges(methanol)

# Print charges
print('METHANOL CHARGES')
print("----------------------------")
print(f'{"Number":<8s}{"Atom":<8s}{"Charge":>12s}')
print("----------------------------")
for i, q in enumerate(methanol_charges):
    element = methanol.elements()[i].symbol()
    charge = q.inUnitsOf(elementary_charge)
    print(f'{i+1:<8d}{element:<8s}{charge:>12.4f}')
print("----------------------------")

# --------------------------------------------------------------------------------------------------
# Create A Box Of Methanol Molecules
# --------------------------------------------------------------------------------------------------
# Define the volume in which the molecules should be packed
u1, u2, u3 = numpy.diag([40.0] * 3) * Angstrom
cell = UnitCell(u1, u2, u3)

# Set the size of the buffer region
buffer_size = 2.0 * Angstrom

# Set Packmol parameters
tolerance = 2.0 * Angstrom
maximum_loops = 30

# Run the Packmol script to generate the packed configuration
methanol_box, packing_successful, packmol_message = runPackmol(
    [(methanol, 900)],
    cell,
    tolerance,
    buffer_size,
    maximum_loops,
)

# --------------------------------------------------------------------------------------------------
# Add A Dreiding Calculator With The Calculated Charges
# --------------------------------------------------------------------------------------------------
# Repeat the charges for each molecule
methanol_box_charges = numpy.tile(
    methanol_charges.inUnitsOf(elementary_charge),
    900
) * elementary_charge

# Create the potential builder
potential_builder = DreidingPotentialBuilder()

# Create the calculator, setting the atomic charges
calculator = potential_builder.createCalculator(
    methanol_box,
    atomic_charges=methanol_box_charges
)

# Assign the calculator to the configuration
methanol_box.setCalculator(calculator)

# --------------------------------------------------------------------------------------------------
# Save The Configuration
# --------------------------------------------------------------------------------------------------
nlsave('Methanol_Box.hdf5', methanol_box)

calculating_charges_example.py

Notes

The ReaxFFAtomicCharges class implements the charge equilibration scheme included in the ReaxFF forcefield [1]. The general charge equilibration method assumes that the change of energy in response to charge polarization can be given as a simple Taylor series expansion, such that the total energy of the molecule, \(E\), can be given as:

\[E(\mathbf{q}) = \sum_i^N \chi_iq_i + \eta_i q_i^2 + \sum_{j<i}^N \frac{q_iq_j}{4\pi \epsilon_0 r_{ij}}\]

Here \(\mathbf{q}\) is the vector of atomic charges, \(\chi_i\) is the atomic electronegativity, \(\eta_i\) is the atomic hardness and \(r_{ij}\) is the inter-atomic distance. Getting the derivative of the energy with respect to an atomic charge gives the chemical potential for that atom. Assuming that the chemical potential of each atom is equal and that the charges sum to a particular value gives a series of linear equations that can be solved for the atomic charges. This gives a simple and yet reasonably effective method for estimating the atomic partial charges.

One of the main modifications of this idea in the ReaxFF charge equilibration scheme is that here the atomic ionization data are treated as adjustable parameters. These are typically tuned to replicate specific charges, and are therefore not necessarily physically meaningful. The adjustable nature of these parameters also means that there are a number of ReaxFF potential sets with parameters tuned to different applications. This is distinct from the QEq charge equilibration method where the ionization data generally retains their atomic values.

Another important modification used in ReaxFF charge equilibration is that the charge-charge electrostatics are calculated with the shielded Coulomb potential

\[V(r_{ij}; \gamma_{ij}) = \frac{1}{4\pi \epsilon_0} \frac{q_iq_j}{\left[r^3+\left(1/\gamma_{ij}\right)^3\right]^{1/3}}\]

This introduces a pairwise shielding parameter \(\gamma_{ij}\). This is calculated as the geometric mean of atomic shielding parameters. The charge-charge interactions are also multiplied by a 7th order polynomial tapering function. This ensures that the interactions are reduced to zero at the cutoff. In periodic system the combination of adjustable atomic parameters, shielded Coulomb function and tapering function allow the charge-charge interactions to converge without Ewald summation. This means that the specified interaction distance cutoff can be used in both molecular and bulk configurations.

To calculate the partial atomic charges for a configuration, an instance of the ReaxFFAtomicCharges class first has to be created. This class has a few parameters to control how it is instantiated. The cutoff used to calculate the charge-charge interactions is controlled with the parameter r_cut. This defines the distance beyond which electrostatic interactions are ignored.

In the ReaxFF scheme each atom has three parameters, the atomic electronegivity, hardness and electrostatic shielding. These can be modified using the parameters chi, eta and shielding respectively. These parameters allow different ReaxFF potential sets to be used. By default the ReaxFFAtomicCharges class uses the potential set of van Duin et. al. [2]. This potential set supports H, C, N, O, Mg, Al, Si, S, K, Cu, Sr and Cs atoms. Each parameter takes dictionary of PeriodicTableElement and PhysicalQuanity of appropriate units. If parameters are not specified for an atom the default atomic parameters are used. The following example shows how to explicitly set the default atomic parameters for the methanol example above.

charger = QEqAtomicCharges(
   chi={Hydrogen: 3.7248*eV, Carbon: 5.9666*eV, Oxygen: 8.5000*eV},
   eta={Hydrogen: 9.6093*eV, Carbon: 7.0000*eV, Oxygen: 8.3122*eV},
   shielding={Hydrogen: 0.8203/Ang, Carbon: 0.9000/Ang, Oxygen: 1.0898/Ang}
)
methanol_charges = charger.charges(methanol)

Once the ReaxFFAtomicCharges class is instantiated, the atomic charges can be calculated using the charges method. This takes the configuration of the system and returns an array of atomic charges. The total charge of the configuration can be set using the total_charge argument. In cases where the charges are sought for a group of individual molecules, the charges can be calculated individually for each molecule using the decompose_bulk argument. This is useful in cases such as a polymer melt or a molecular liquid. Calculating the charges of each molecule individually ensures that each molecule has the correct total charge, as well as removing any artifacts of how the molecules are packed. It also can speed up the calculation of the charges significantly. This is currently only supported for BulkConfiguration configurations.