QEqAtomicCharges

class QEqAtomicCharges(r_cut=PhysicalQuantity(10.0, Ang), error=PhysicalQuantity(1e-05, 1 / Ang), chi=None, eta=None, shielding=None, max_cycles=50, mixing=0.5)

Enables the calculation of atomic partial charges using the QEq charge equilibration scheme of Rappe and Goddard.

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

  • error (PhysicalQuantity of type inverse length) – The error in the 1/r lattice sum.
    Default: 0.00001 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 QEq method.

  • 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 QEq method.

  • 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 QEq method.

  • max_cycles (int) – The maximum number of iterations performed in QEq calculations.
    Default: 50

  • mixing (float) – The linear mixing between new and old charges calculated by QEq.
    Default: 0.5

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 QEq 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 = QEqAtomicCharges()
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 QEqAtomicCharges class implements the QEq charge equilibration scheme of Rappe and Goddard [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 modification of this idea in the QEq scheme is that the charge-charge interaction term in the above equation is not treated with point charges, but rather as the electrostatic interactions of two normalized Slater functions. This leads to a better approximation of charge shielding at short ranges. The Slater functions on each atom can be given as:

\[\phi(\mathbf{r}; n,\zeta) = Nr^{n-1}e^{-\zeta r}\]

Here \(N\) is a normalization constant, \(\zeta\) is the principal quantum number and \(\zeta\) is the orbital exponent. This sets the effective size of the Slater function on the atom. As Slater electrostatic integrals are difficult to calculate, a number of effective approximations of this term have been suggested [2].

Periodic systems also present an additional challenge to the QEq scheme. Due to the long-range nature of electrostatic interactions, in periodic systems the charge-charge interactions should be calculated with Ewald summation [3]. As the Slater electrostatic term approaches \(\frac{q_iq_j}{4\pi \epsilon_0 r_{ij}}\) at long distances, this long range interaction can be replaced with the appropriate Ewald sum.

Another important modification in QEq is in how hydrogen is treated. Since the atomic orbitals of hydrogen are significantly different to its molecular orbitals, it is not sufficient to use atomic parameters in calculating charge equilibration in molecules. Charge dependent corrections to the atomic hardness and Slater orbital exponent are therefore used to make hydrogen respond to charge polarization in a more molecular way. As these corrections are charge dependent this means that the QEq charge calculations have to be solved self consistently. This can be done using a simple linear mixing scheme.

To calculate the partial atomic charges for a configuration, an instance of the QEqAtomicCharges 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. In the case of a MoleculeConfiguration this defines the distance beyond which electrostatic interactions are ignored. For large molecules such as polymers this may need to be increased so that all interactions are considered. In the case of a BulkConfiguration this defines the distance beyond which electrostatic interactions are assumed to be \(\frac{q_iq_j}{4\pi \epsilon_0 r_{ij}}\) and calculated with Ewald summation. Before this distance electrostatic interactions are assumed to be Slater electrostatic integrals.

In the QEq scheme each atom has three parameters, the atomic electronegativity, hardness and Slater orbital exponent. These can be modified using the parameters chi, eta and shielding respectively. This can be important in cases where atomic charges are being sought for ionic systems. In the QEq scheme the energy of the atom is expanded around a neutral charge. In ionic systems it is more appropriate to expand the energy around the ionic charge of the atom. This different atomic reference state is reflected in different atomic parameters. 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: 4.528*eV, Carbon: 5.343*eV, Oxygen: 8.741*eV},
   eta={Hydrogen: 6.9452*eV, Carbon: 5.063*eV, Oxygen: 6.682*eV},
   shielding={Hydrogen: 2.0216/Ang, Carbon: 1.6469/Ang, Oxygen: 1.8685/Ang}
)
methanol_charges = charger.charges(methanol)

Finally there are two parameters that control the iterative solving of the QEq equations. This is only necessary in configurations that contain hydrogen. The mixing parameter takes a number between 0 and 1 that represents the amount of charges to add from the current and previous solutions to add at each iteration. The parameter max_cycles sets the maximum number of self consistency iterations. In some cases these parameters may need to be modified to help the self consistent solution to converge.

Once the QEqAtomicCharges 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.