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)
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:
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:
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.