Skip to content

eeq

build_a_matrix(atomic_numbers, positions)

Build the A matrix for EEQ charge calculation.

Parameters:

Name Type Description Default
atomic_numbers ndarray

Array of atomic numbers for each atom

required
positions ndarray

Array of atomic positions (shape: N x 3)

required

Returns:

Type Description

np.ndarray: A matrix for the EEQ calculation

Source code in chmpy/core/eeq.py
def build_a_matrix(atomic_numbers, positions):
    """
    Build the A matrix for EEQ charge calculation.

    Args:
        atomic_numbers (np.ndarray): Array of atomic numbers for each atom
        positions (np.ndarray): Array of atomic positions (shape: N x 3)

    Returns:
        np.ndarray: A matrix for the EEQ calculation
    """
    N = len(atomic_numbers)
    sqrt_pi_fac = np.sqrt(2.0 / np.pi)
    A = np.zeros((N + 1, N + 1))

    positions_bohr = positions * ANGSTROM_TO_BOHR

    r = pdist(positions_bohr)
    r2 = r**2

    widths = WIDTH[atomic_numbers]

    i, j = np.triu_indices(N, k=1)

    ri_squared = widths[i] ** 2
    rj_squared = widths[j] ** 2
    gamma = 1.0 / (ri_squared + rj_squared)

    values = erf(np.sqrt(r2 * gamma)) / np.sqrt(r2)

    A_temp = np.zeros((N, N))
    A_temp[i, j] = values
    A_temp[j, i] = values  # Mirror values (symmetric matrix)

    A[:N, :N] = A_temp

    diagonal_values = ETA[atomic_numbers] + sqrt_pi_fac / WIDTH[atomic_numbers]
    np.fill_diagonal(A[:N, :N], diagonal_values)

    A[N, :N] = 1.0
    A[:N, N] = 1.0
    A[N, N] = 0.0

    return A

build_x_vector(atomic_numbers, cn, charge=0.0)

Build the X vector for EEQ charge calculation.

Parameters:

Name Type Description Default
atomic_numbers ndarray

Array of atomic numbers for each atom

required
cn ndarray

Array of coordination numbers for each atom

required
charge float

Total charge of the system

0.0

Returns:

Type Description

np.ndarray: X vector for the EEQ calculation

Source code in chmpy/core/eeq.py
def build_x_vector(atomic_numbers, cn, charge=0.0):
    """
    Build the X vector for EEQ charge calculation.

    Args:
        atomic_numbers (np.ndarray): Array of atomic numbers for each atom
        cn (np.ndarray): Array of coordination numbers for each atom
        charge (float): Total charge of the system

    Returns:
        np.ndarray: X vector for the EEQ calculation
    """
    N = atomic_numbers.shape[0]
    eps = 1e-14  # Avoid singularity with 0
    X = np.empty(N + 1)
    X[:N] = -CHI[atomic_numbers] + cn * KCN_PARAM[atomic_numbers] / np.sqrt(cn + eps)
    X[N] = charge
    return X

calculate_coordination_numbers(atomic_numbers, positions)

Calculate coordination numbers for all atoms in a molecule or crystal.

Parameters:

Name Type Description Default
atomic_numbers ndarray

Array of atomic numbers for each atom

required
positions ndarray

Array of atomic positions (shape: N x 3)

required

Returns:

Type Description

np.ndarray: Array of coordination numbers for each atom

Source code in chmpy/core/eeq.py
def calculate_coordination_numbers(atomic_numbers, positions):
    """
    Calculate coordination numbers for all atoms in a molecule or crystal.

    Args:
        atomic_numbers (np.ndarray): Array of atomic numbers for each atom
        positions (np.ndarray): Array of atomic positions (shape: N x 3)

    Returns:
        np.ndarray: Array of coordination numbers for each atom
    """
    N = len(atomic_numbers)
    kcn_value = 7.5  # Constant from the C++ implementation

    cn = np.zeros(N)
    cutoff = 25.0  # Cutoff distance in Angstroms squared

    dists = pdist(positions) * ANGSTROM_TO_BOHR

    i, j = np.triu_indices(N, k=1)

    mask = dists <= cutoff
    dists = dists[mask]
    i_filt, j_filt = i[mask], j[mask]

    rc = (
        COVALENT_D3[atomic_numbers[i_filt]] + COVALENT_D3[atomic_numbers[j_filt]]
    ) * ANGSTROM_TO_BOHR

    counts = 0.5 * (1.0 + erf(-kcn_value * (dists / rc - 1.0)))

    np.add.at(cn, i_filt, counts)
    np.add.at(cn, j_filt, counts)

    return cn

calculate_eeq_charges(atomic_numbers, positions, charge=0.0)

Calculate EEQ partial charges for a set of atoms.

Parameters:

Name Type Description Default
atomic_numbers ndarray

Array of atomic numbers for each atom

required
positions ndarray

Array of atomic positions (shape: N x 3)

required
charge float

Total charge of the system

0.0

Returns:

Type Description

np.ndarray: Array of partial charges for each atom

Source code in chmpy/core/eeq.py
def calculate_eeq_charges(atomic_numbers, positions, charge=0.0):
    """
    Calculate EEQ partial charges for a set of atoms.

    Args:
        atomic_numbers (np.ndarray): Array of atomic numbers for each atom
        positions (np.ndarray): Array of atomic positions (shape: N x 3)
        charge (float): Total charge of the system

    Returns:
        np.ndarray: Array of partial charges for each atom
    """
    # Calculate coordination numbers
    cn = calculate_coordination_numbers(atomic_numbers, positions)
    A = build_a_matrix(atomic_numbers, positions)
    X = build_x_vector(atomic_numbers, cn, charge)

    Q = np.linalg.solve(A, X)
    return Q[:-1]