RepeatESPPartialCharges

class RepeatESPPartialCharges(configuration, potential, total_charge=None, tag_prefix=None, charge_weight=None, inner_vdw_shell=None, outer_vdw_shell=None, stride=None, r_cut=None, damping_factor=None)

Calculate atomic partial charges to reproduce the electrostatic potential based on the REPEAT algorithm.

Parameters:

• configuration (MoleculConfiguration | BulkConfiguration) – The configuration for which the charges are fitted.

• potential (ElectrostaticDifferencePotential) – The electrostatic potential the charges are fitted to.

• total_charge (PhysicalQuantity of type charge) – The total charge of the configuration.

• tag_prefix (str) – Prefix of the tags on the configuration used to constrain charges.

• charge_weight (float | list of float) – Weights used to constrain fitted charges to the partial charges on the configuration.

• inner_vdw_shell (float) – Ratio of the van der Waals radius of the atoms used to exclude grid points close to an atom.

• outer_vdw_shell (float) – Ratio of the van der Waals radius of the atoms used to exclude grid points far away from any atom.

• stride (int) – Stride used to skip grid points in the electrostatic potential. A value of 1 specifies that every grid point is used.

• r_cut (PhysicalQuantity of type length) – The cutoff used in calculating electrostatic interactions in BulkConfiguration.

• damping_factor (PhysicalQuantity of type inverse length) – The damping factor used in calculating electrostatic interactions in BulkConfiguration.

dampingFactor()

Returns:

Damping factor used for fitted electrostatic interactions in bulk materials.

Return type:

PhysicalQuantity of type inverse length

evaluate(atom_index=None)

Return the partial charges of the selected atom or all atoms.

Parameters:

atom_index (int) – Index of the selected atom.
Default: All atoms.

Returns:

A single charge if atom_index is not None, otherwise an array of charges.

Return type:

PhysicalQuantity of type charge

groupIndices()

Returns:

Indices of groups of atoms constrained to have the same fitted charge.

Return type:

list of lists

innerVdwShell()

Returns:

The inner radius factor of the shell for electrostatic potential grid points.

Return type:

float

metatext()

Returns:

The metatext of the object or None if no metatext is present.

Return type:

str | None

nlprint(stream=None)

Print a string containing an ASCII table useful for plotting the AnalysisSpin object.

Parameters:

stream (python stream) – The stream the table should be written to.
Default: NLPrintLogger()

outerVdwShell()

Returns:

The outer radius of the shell for electrostatic potential grid points.

Return type:

float

partialCharges()

Returns:

The fitted partial charges for the molecule.

Return type:

PhysicalQuantity of type charge

rCut()

Returns:

Electrostatic potential distance cutoff used in bulk materials.

Return type:

PhysicalQuantity of type length

relativeRootMeanSquareError()

Returns:

The relative root mean squared error of the electrostatic potential fit.

Return type:

float

restrainCharges()

Returns:

The reference charges used to weight the fitted charges.

Return type:

PhysicalQuantity of type charge

restrainWeights()

Returns:

The weights used to restrain the fitted charges.

Return type:

ndarray

setMetatext(metatext)

Set a given metatext string on the object.

Parameters:

metatext (str | None) – The metatext string that should be set. A value of “None” can be given to remove the current metatext.

stride()

Returns:

Density of electrostatic potential points selected for fitting partial charges.

Return type:

int

tags()

Returns:

Tags used to constrain atoms to have the same partial charges.

Return type:

list of str

totalCharge()

Returns:

The total charge of the configuration the charges are being fitted to.

Return type:

PhysicalQuantity of type charge

uniqueString()

Return a unique string representing the state of the object.

Usage Examples

Assign partial charges to the battery electrolyte molecule ethylene carbonate. In this example the OPLS potential is used to constrain atoms with the same type to have the same partial charges.

# -*- coding: utf-8 -*-
setVerbosity(MinimalLog)

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

# Define coordinates
cartesian_coordinates = [[ 1.609431413032,  0.050029424266, -0.200214465477],
                         [ 2.821872471776,  0.074345609581, -0.33257452873 ],
                         [ 0.786730069642,  0.285844748805, -1.270661995358],
                         [ 1.046027032925, -0.213547522321,  1.021790323208],
                         [-0.320536198865, -0.108625807588,  0.763421170161],
                         [-0.484549258604,  0.072564433479, -0.73758247059 ],
                         [-0.839575976013, -1.027933189299,  1.111755833993],
                         [-0.732343636487,  0.769867954492,  1.307351219363],
                         [-0.918484269576, -0.844342167145, -1.192644130163],
                         [-1.139565310547,  0.94179651573 , -0.964104219255]]*Angstrom

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

# Add tags
molecule_configuration.addTags('OPLS_CT_226', [4, 5])
molecule_configuration.addTags('OPLS_C_32',   [0])
molecule_configuration.addTags('OPLS_HC_63',  [6, 7, 8, 9])
molecule_configuration.addTags('OPLS_OS_22',  [2, 3])
molecule_configuration.addTags('OPLS_O_23',   [1])

# Add bonds
bonds = [[0, 1],
         [0, 2],
         [0, 3],
         [2, 5],
         [3, 4],
         [4, 5],
         [4, 6],
         [4, 7],
         [5, 8],
         [5, 9]]
molecule_configuration.setBonds(bonds)

# -------------------------------------------------------------
# Calculator
# -------------------------------------------------------------
numerical_accuracy_parameters = NumericalAccuracyParameters(
    density_mesh_cutoff=125.0*Hartree,
)

calculator = LCAOCalculator(
    numerical_accuracy_parameters=numerical_accuracy_parameters,
)

molecule_configuration.setCalculator(calculator)
nlprint(molecule_configuration)
molecule_configuration.update()
nlsave('ethylene_earbonate_charges.hdf5', molecule_configuration)

# -------------------------------------------------------------
# Optimize Geometry
# -------------------------------------------------------------
molecule_configuration = OptimizeGeometry(
    molecule_configuration,
    max_forces=0.05*eV/Ang,
    max_steps=200,
    max_step_length=0.2*Ang,
    trajectory_filename='ethylene_earbonate_charges_trajectory.hdf5',
    trajectory_interval=5.0*Minute,
    restart_strategy=RestartFromTrajectory(),
    optimize_cell=False,
    optimizer_method=LBFGS(),
    enable_optimization_stop_file=True,
)
nlsave('ethylene_earbonate_charges.hdf5', molecule_configuration)
nlprint(molecule_configuration)

# -------------------------------------------------------------
# Electrostatic Difference Potential
# -------------------------------------------------------------
electrostatic_difference_potential = ElectrostaticDifferencePotential(molecule_configuration)
nlsave('ethylene_earbonate_charges.hdf5', electrostatic_difference_potential)

# -------------------------------------------------------------
# Repeat ESP Partial Charges
# -------------------------------------------------------------
repeat_esp_partial_charges = RepeatESPPartialCharges(
    molecule_configuration,
    electrostatic_difference_potential,
    total_charge=0*elementary_charge,
    inner_vdw_shell=1.0,
    outer_vdw_shell=3.0,
    stride=1,
    tag_prefix='OPLS_',
    )
molecule_configuration.setPartialCharges(repeat_esp_partial_charges.evaluate())
nlsave('ethylene_earbonate_charges.hdf5', molecule_configuration)
nlsave('ethylene_earbonate_charges.hdf5', repeat_esp_partial_charges)

ethylene_carbonate_charges.py

Notes

The RespeatESPPartialCharges class fits atomic partial charges to either a MoleculeConfiguration or BulkConfiguration based on a given electrostatic potential. This potential is typically calculated with density functional theory (DFT). Theseo partial charges can then be used in a forcefield based calculation to model electrostatic interactions between atoms. As the charges are fit to a DFT electrostatic potential, it allows the long-range electrostatic interactions to be similar to those predicted by DFT, resulting in a more accurate simulation.

The fitting processes uses the REPEAT algorithm to assign the partial charges[1]. This fits the charges using a linear least squares approach, minimizing the difference between the supplied electrostatic potential and the potential calculated from the fitted partial charges. The grid points in the electrostatic potential used in the fitting are selected such that they are all within a shell around the molecule. Grid points that are too close to an atom are within the electron density around the atom where the simple approximation of point charges breaks down. Similarly the electrostatic potential far from any atom is weakly influenced by many surrounding atoms and therefore not useful in isolating the contribution of single atoms to the electrostatic potential. The inner and outer radii used to define the shell of grid points used in the fitting can be set using the options inner_vdw_shell and outer_vdw_shell. These are scale factors that are combined with the van der Waals radius of each atom to determine the shell radii. Depending on the density mesh cutoff used in the DFT calculation, the electrostatic potential can also have a large number of closely spaced grid points. It is possible to down-sample the electrostatic potential grid, skipping over a set number of points in each direction. This is done using the stride argument, which takes an integer describing the increment used when reading points from the grid. The default value of 1 uses all the grid points in the given electrostatic potential.

As the charges are fit using a simple least squares method, it is possible that some atoms may be assigned chemically unreasonable charges, such as charges exceeding the number of valence electrons. This is most commonly an issue in configurations that have buried atoms, that is atoms that are not near any electrostatic grid points used in the fitting. This may happen for instance in atoms contained within a bulk surface. To help charges retain chemically meaningful values, some constraints can be added to the fit. One constraint that can be added is to constrain groups of atoms to have the same partial charge. This is done using tags. Tags identifying atom that have the same charge can be added to the configuration with a set prefix, with one tag for each atom. The tag prefix is then given as input to the RespeatESPPartialCharges class using the tag_prefix argument. Tags from bonded forcefields such as DREIDING or OPLS can be used, or custom tags can also be set. The given example with ethylene carbonate uses the OPLS tags to constrain atoms with the same OPLS type to have the same partial charge.

It is also possible to set reference charges for the fit, and weight the fitting towards these charges. This weighting can give reasonable values for charges that are not well determined by the electrostatic potential, while still allowing some influence from the potential. Reference charges are taken from the existing atomic partial charges on the configuration. Weights for each atom are given with the argument charge_weight. If a single value is given, this value is used for all atoms. Giving individual values for each atom allows only under determined charges to be weighted to particular values. The size of the weight determines how strongly the existing partial charge influences the values of the fitted charges.

In molecular and bulk configurations the summation of the fitted electrostatic potential is done in slightly different ways. In molecules the fitted electrostatic potential is represented with the normal Coulomb potential, which does not require additional parameters. In the case of bulk configurations the fitted electrostatic potential is represented with a damped shifted force potential[2]. This requires both a distance cutoff and a damping factor. These can be modified using the arguments r_cut and damping_factor respectively.

[1]

Carlos Campana, Bastien Mussard, and Tom K. Woo. Electrostatic potential derived atomic charges for periodic systems Using a modified error functional. Journal of Chemical Theory and Computation, 5(10):2866–2878, 2009. URL: https://pubs.acs.org/doi/10.1021/ct9003405, doi:https://doi.org/10.1021/ct9003405.

[2]

De-Li Chen, Abraham C. Stern, Brian Space, and J. Karl Johnson. Atomic charges derived from electrostatic potentials for molecular and periodic systems. Journal of Physical Chemistry A, 114(37):10225–10233, 2009. URL: https://pubs.acs.org/doi/10.1021/jp103944q, doi:https://doi.org/10.1021/jp103944q.