pyef.analysis

Electrostatics Analysis Module for Quantum Chemistry Calculations.

This module provides the Electrostatics class for computing electrostatic properties (ESP, electric fields, partial charges) from quantum chemistry calculations. It interfaces with Multiwfn for charge analysis and supports various charge partitioning schemes (Hirshfeld, Hirshfeld-I, CHELPG, etc.).

Key Features:

  • Multiple charge partitioning schemes

  • ESP and electric field calculations

  • Multipole moment analysis

  • QMMM point charge integration

  • Dielectric environment effects

  • ECP (Effective Core Potential) support

Dependencies:

  • Multiwfn (external program)

  • NumPy, pandas, scipy

  • OpenBabel, BioPandas

Module Contents

Classes

Electrostatics

Compute electrostatic properties from quantum chemistry calculations.

Functions

visualize_charges_pdb(xyz_file, charges, pdb_name)

Create a PDB file with charges visualized in the b-factor column.
class Electrostatics(molden_paths=None, xyz_paths=None, esp_atom_idx=None, ptchg_paths=None, pdb_charge_paths=None, **kwargs)

Compute electrostatic properties from quantum chemistry calculations.

This class processes series of quantum chemistry calculations to extract electrostatic properties including ESP, electric fields, and partial charges using various charge partitioning schemes via Multiwfn.

molden_paths

Paths to .molden files, one per job (absolute or relative).

Type:

list of str

xyz_paths

Paths to .xyz files, one per job (absolute or relative).

Type:

list of str

esp_atom_idx

Atom indices for ESP evaluation (typically metal centers).

Type:

list of int

config
Configuration dictionary containing:

  • hasECP (bool): Whether ECP was used in calculation

  • includePtChgs (bool): Include point charges in ESP calculation

  • ptChgfp (str): Path to point charge file

  • molden_filename (str): Name of molden file (compatibility)

  • xyzfilename (str): Name of xyz coordinate file (compatibility)

  • rerun (bool): Force recalculation of charges

  • maxIHirshBasis (int): Max basis functions for Hirshfeld-I

  • maxIHirshFuzzyBasis (int): Max basis for fuzzy Hirshfeld-I

  • ECP (str): ECP basis set family

  • dielectric (float): Dielectric constant

  • dielectric_scale (float): Scaling factor for dielectric

  • excludeAtomfromEcalc (list): Atoms to exclude from E-field calc

  • changeDielectBoundBool (bool): Use dielectric=1 for bonded atoms

  • visualize_ef (bool): Create PDB files for atom-wise E-fields

  • visualize_charges (bool): Create PDB files for partial charges

  • visualize_per_bond (bool): Create separate PDB files for each bond

  • skip_missing_files (bool): Skip calculations if required files missing (default: False)

Type:

dict

dict_of_calcs

Mapping of charge scheme names to Multiwfn command codes.

Type:

dict

dict_of_multipole

Mapping of multipole schemes to Multiwfn command sequences.

Type:

dict

ptchgdf

DataFrame containing point charge data for QMMM calculations.

Type:

pd.DataFrame or None

periodic_table

Mapping of element symbols to atomic numbers.

Type:

dict

amassdict

Atomic properties (mass, atomic number, covalent radius, valence).

Type:

dict

Examples

>>> molden_paths = ['calc1/scr/optim.molden', 'calc2/scr/optim.molden']
>>> xyz_paths = ['calc1/scr/optim.xyz', 'calc2/scr/optim.xyz']
>>> es = Electrostatics(molden_paths, xyz_paths, dielectric=2.0)
>>>
>>> # For ESP calculations (metal indices required)
>>> es_with_metals = Electrostatics(molden_paths, xyz_paths, esp_atom_idx=[0, 0], dielectric=2.0)
>>> data = es_with_metals.getESP(['Hirshfeld', 'Hirshfeld_I'], 'esp_results', ...)
run_multiwfn(command, input_commands, output_file=None, description='Multiwfn calculation', capture_output=True)

Centralized Multiwfn runner with proper error handling and output display.

Parameters:

  • command (str) – The shell command to run Multiwfn (e.g., “multiwfn file.molden”)

  • input_commands (list of str) – List of commands to send to Multiwfn’s stdin

  • output_file (str, optional) – Path to redirect stdout to a file

  • description (str, optional) – Description of the calculation for error messages

  • capture_output (bool, optional) – If True, capture and display output. If False, let it stream to terminal

Returns:

(stdout, stderr, returncode) from the Multiwfn process

Return type:

tuple

Raises:

RuntimeError – If Multiwfn fails with non-zero exit code

updateCalcSettings(key, value)

Update calculation settings dictionary.

Parameters:

  • key (str) – Setting key to update.

  • value (any) – New value for the setting.

setExcludeAtomFromCalc(lstExcludeAtoms)

Set atoms to exclude from electric field calculations.

Parameters:

lstExcludeAtoms (list of int) – Atom indices (0-indexed) to exclude from E-field calculations.

includePtChgs(name_ptch_file)

Include point charges in ESP calculation (for QMMM).

Parameters:

name_ptch_file (str or list) – If str, uses a single filename (e.g., ‘pointcharges.txt’) that is searched for in each job directory. If list, uses per-job file paths, one per job; use None for jobs without point charges.

Notes

Sets includePtChgs config to True and stores filename/paths. Point charges are typically from MM region in QM/MM calculations.

For filename mode (string), the file is automatically located in each job’s working directory. For list mode, paths may be absolute or relative to the current working directory.

Examples

# Single filename (applies per job directory) estat.includePtChgs(‘pointcharges.txt’)

# Per-job paths estat.includePtChgs([‘/path/job1/ptchg.txt’, None, ‘/path/job3/ptchg.txt’])

set_dielec_scale(dielec)

Set dielectric scaling factor for point charges.

Parameters:

dielec (float) – Scaling factor for dielectric in QMMM calculations.

initialize_excludeAtomsFromEfieldCalc(atom_to_exclude)

Initialize atoms to exclude from electric field calculation.

Parameters:

atom_to_exclude (list of int) – List of atom indices (0-indexed) to exclude.

minDielecBonds(bool_bonds)

Set dielectric to 1 for bonded atoms in ESP calculation.

Parameters:

bool_bonds (bool) – If True, use dielectric=1 for atoms bonded to ESP center. Helps avoid over-estimating screening from bound atoms.

runlowmemory()

Enable low-memory mode for large systems (>300 atoms).

Notes

Switches Hirshfeld-I calculation to slower but lower-memory algorithm. Useful for large systems that would otherwise exceed memory limits.

changeDielectric(dlc)

Change dielectric constant of the environment.

Parameters:

dlc (float) – Dielectric constant of solvent/environment.

Notes

Used to model solvation effects in ESP calculations. Default is 1 (vacuum). Common values: water=80, protein=4.

set_molden_filename(new_name)

Set the expected molden filename.

Parameters:

new_name (str) – Name of the molden file to read.

set_xyzfilename(new_name)

Set the expected XYZ coordinate filename.

Parameters:

new_name (str) – Name of the XYZ file to read.

rePrep()

Re-run data preparation (file extraction and validation).

Notes

Useful if files have been updated or regenerated.

fix_allECPmolden()

Fix molden files that used Effective Core Potentials.

Notes

ECPs can cause artifacts in molden files that make them incompatible with Multiwfn. This method reformats the molden files to fix these issues. Processes all molden/xyz file pairs.

prepData()

Validate that all specified .molden and .xyz files exist.

Checks that all file paths provided during initialization are valid and accessible. Removes invalid entries from processing lists.

Raises:

FileNotFoundError – If any of the specified .molden or .xyz files cannot be found.

prepData_pdb()

Validate PDB charge files for pdb_only mode.

Checks that all PDB files provided during initialization exist. Removes invalid entries from processing lists.

Raises:

FileNotFoundError – If no valid PDB charge files are found.

static getmultipoles(multipole_name)

Parse multipole moments from Multiwfn output file.

Parameters:

multipole_name (str) – Path to multipole moment output file from Multiwfn.

Returns:

List of dictionaries containing multipole data for each atom: - ‘Index’: atom index (str) - ‘Element’: element symbol (str) - ‘Atom_Charge’: atomic charge (float) - ‘Dipole_Moment’: dipole moment vector [x, y, z] (list of float) - ‘Quadrupole_Moment’: 3x3 quadrupole tensor (np.ndarray)

Return type:

list of dict

Notes

Parses Multiwfn output format with sections separated by asterisks. Extracts monopole (charge), dipole, and quadrupole moments for each atom.

getPtChgs(filename_pt)

Load point charges from file (for QMMM calculations).

Parameters:

filename_pt (str) – Path to point charge file (typically from Amber/CHARMM).

Returns:

DataFrame with columns: [‘Atom’, ‘charge’, ‘x’, ‘y’, ‘z’] All atoms labeled as ‘pnt’ (point charge).

Return type:

pd.DataFrame

Notes

File format: whitespace-delimited with 2 header lines to skip.

getPdbCharges(pdb_path)

Load point charges from a PDB file using the B-factor column.

Parameters:

pdb_path (str) – Path to PDB file containing charges in the B-factor column.

Returns:

DataFrame with columns: [‘Atom’, ‘charge’, ‘x’, ‘y’, ‘z’, ‘ResName’]. Row position equals the 0-indexed atom index used throughout pyEF (e.g. in input_bond_indices, exclude_atoms, esp_atom_idx).

Return type:

pd.DataFrame

Notes

Reads ATOM/HETATM records and parses fixed-width PDB columns. The B-factor column is interpreted as the partial charge (e). ‘ResName’ is the standard PDB residue name field (columns 18-20) and is used by getEfield_byResidue()/getESP_byResidue() to group atoms.

static mapcount(filename)

Rapidly count lines in a file using memory mapping.

Parameters:

filename (str) – Path to file to count.

Returns:

Number of lines in the file.

Return type:

int

Notes

Uses mmap for efficient line counting in large files.

monopole_esp(path_to_xyz, espatom_idx, charge_range, charge_file, ptchg_path=None)

Calculate electrostatic potential at specified atom.

Parameters:

  • path_to_xyz (str) – Path to XYZ coordinate file.

  • espatom_idx (int) – Atom index (0-indexed) where ESP will be calculated.

  • charge_range (list of int) – Indices of atoms to include in ESP calculation.

  • charge_file (str) – Path to partial charge file from Multiwfn.

  • ptchg_path (str or None, optional) – Path to point charge file. If None, no point charges are included.

Returns:

[ESP (float), atom_symbol (str)] ESP in Volts at the specified atom.

Return type:

list

Notes

  • Uses Coulomb’s law: V = k * q / r with dielectric screening

  • Point charges from QMMM included if ptchg_path is provided

  • Bonded atoms treated with dielectric=1 if config[‘changeDielectBoundBool’] is True

  • Units: Volts (N·m/C = J/C)

monopole_esp_df(df_charges, espatom_idx, charge_range)

Calculate ESP at a specified atom using charges from a DataFrame.

This is the pdb_only equivalent of monopole_esp(): it reads directly from a DataFrame (e.g. produced by getPdbCharges) instead of a file. The changeDielectBoundBool feature is not supported in this mode because bond detection requires an XYZ file.

Parameters:

  • df_charges (pd.DataFrame) – DataFrame with columns [‘Atom’, ‘x’, ‘y’, ‘z’, ‘charge’]. Typically returned by getPdbCharges().

  • espatom_idx (int) – Atom index (0-indexed) at which the ESP is calculated.

  • charge_range (list of int) – Indices of atoms to include in the ESP sum.

Returns:

[ESP (float in Volts), atom_symbol (str)]

Return type:

list

calc_firstTermE(charge_range, charge_file)

Input: espatom_idx: integer of atom index charge_range: list of integers of atom indices charge_file: string of charge filename Output: list of E-field vector and atomic symbol

monopole_atomicEfield(espatom_idx, charge_range, charge_file, df_ptchg=None)

Input: espatom_idx: integer of atom index charge_range: list of integers of atom indices charge_file: string of charge filename df_ptchg: optional pre-loaded point charge dataframe Output: list of E-field vector and atomic symbol

pointcharge_atomicEfield(espatom_idx, charge_range, df_charges)

Calculate E-field using point charges from a DataFrame only.

Parameters:

  • espatom_idx (int) – Atom index (0-indexed) where the E-field is calculated.

  • charge_range (list or range) – Atom indices to include in the E-field calculation.

  • df_charges (pd.DataFrame) – DataFrame with columns: [‘Atom’, ‘x’, ‘y’, ‘z’, ‘charge’].

Returns:

[E_vec, position_vec, atom_symbol, atom_wise_contributions] E_vec is in V/Angstrom.

Return type:

list

pointcharge_atomicESP(espatom_idx, charge_range, df_charges)

Calculate ESP at an atom using point charges, with atom-wise contributions.

ESP equivalent of pointcharge_atomicEfield(): same total value as monopole_esp_df(), but additionally returns the contribution of each individual atom so they can be summed per residue group.

Parameters:

  • espatom_idx (int) – Atom index (0-indexed) at which the ESP is calculated.

  • charge_range (list or range) – Atom indices to include in the ESP sum.

  • df_charges (pd.DataFrame) – DataFrame with columns: [‘Atom’, ‘x’, ‘y’, ‘z’, ‘charge’].

Returns:

[ESP (float, Volts), atom_symbol (str), atomwise_contributions (np.ndarray, Volts)]

Return type:

list

monopole_pointEfield(point, charge_range, charge_file, df_ptchg=None)

Compute E-field at an arbitrary Cartesian point using monopole charges.

Parameters:

  • point (array-like of float) – [x, y, z] in Angstrom of the evaluation point.

  • charge_range (list or range) – Atom indices to include in the calculation.

  • charge_file (str) – Path to the charge file.

  • df_ptchg (pd.DataFrame or None) – Optional pre-loaded point charge dataframe.

Returns:

[E_vec (V/Angstrom), position_vec (m), atom_wise_contributions (V/Angstrom)]

Return type:

list

pointcharge_pointEfield(point, charge_range, df_charges)

Compute E-field at an arbitrary Cartesian point using point charges from a DataFrame.

Parameters:

  • point (array-like of float) – [x, y, z] in Angstrom of the evaluation point.

  • charge_range (list or range) – Atom indices to include in the calculation.

  • df_charges (pd.DataFrame) – DataFrame with columns: [‘Atom’, ‘x’, ‘y’, ‘z’, ‘charge’].

Returns:

[E_vec (V/Angstrom), position_vec (m), atom_wise_contributions (V/Angstrom)]

Return type:

list

compute_esp(q_vals, probe_coords)

Compute ESP at probe_coords due to point charges at q_coords.

compute_esp_from_qm(qm_charges, mm_coords)

Compute the ESP at MM atoms due to QM charges.

multipole_atomicEfield(idx_atom, charge_range, xyz_file, atom_multipole_file, df_ptchg=None)

Input: idx_atom: integer of atom index charge_range: list of integers of atom indices xyz_file: string of xyz filename atom_multipole_file: string of multipole filename df_ptchg: optional pre-loaded point charge dataframe Output: list of E-field vector and atomic symbol, Efields are reported in units of Volts/Angstrom

multipole_pointEfield(point, charge_range, xyz_file, atom_multipole_file, df_ptchg=None)

Compute E-field at an arbitrary Cartesian point using multipole expansion.

Parameters:

  • point (array-like of float) – [x, y, z] in Angstrom of the evaluation point.

  • charge_range (list or range) – Atom indices to include in the calculation.

  • xyz_file (str) – Path to the XYZ structure file.

  • atom_multipole_file (str) – Path to the multipole file.

  • df_ptchg (pd.DataFrame or None) – Optional pre-loaded point charge dataframe.

Returns:

[E_vec (V/Angstrom), position_vec (m), atom_wise_contributions (V/Angstrom)]

Return type:

list

multipole_esp(xyzfilepath, atom_multipole_file, charge_range, idx_atom, ptchg_path=None)

Input: idx_atom: integer of atom index charge_range: list of integers of atom indices xyfilepath: string of xyz filename atom_multipole_file: string of multipole filename ptchg_path: optional path to point charge file Output: list of ESP and atomic symbol Units of ESP is Volts!

calc_bond_dipoles(bond_indices, xyz_filepath, atom_multipole_file, bool_multipole, charge_df=None)

Calculate bond dipole moments for given bond indices.

For each bond between atoms A and B, calculates the bond dipole moment vector relative to the geometric center (midpoint) of the bond. This is origin-independent for both neutral and charged bonds: μ_bond = q_A * (r_A - r_center) + q_B * (r_B - r_center) which simplifies to: μ_bond = (q_A - q_B) * (r_A - r_B) / 2

Parameters:

  • bond_indices (list of tuples) – List of (atomA_idx, atomB_idx) tuples specifying bonds

  • xyz_filepath (str) – Path to XYZ structure file

  • atom_multipole_file (str) – Path to multipole/charge file

  • bool_multipole (bool) – If True, use multipole file; if False, use monopole charges

  • charge_df (pd.DataFrame or None) – Optional DataFrame with monopole charges (columns: Atom, x, y, z, charge). If provided, overrides charge file input.

Returns:

List of bond dipole information dictionaries containing: - ‘bond_dipole_vec’: [x, y, z] dipole vector in Debye - ‘bond_dipole_mag’: magnitude of dipole in Debye - ‘bond_indices’: (atomA_idx, atomB_idx) - ‘charges’: (charge_A, charge_B)

Return type:

list

bondEfield(bond_indices, xyz_filepath, atom_multipole_file, all_lines, bool_multipole, df_ptchg=None)

Input: bond_indices: list of tuples of atom indices xyz_filepath: string of xyz filename atom_multipole_file: string of multipole filename all_lines: list of strings of lines in xyz file df_ptchg: optional pre-loaded point charge dataframe Output: list of E-field vector and atomic symbol in units of V/Angstrom

bondEfield_pointcharges(bond_indices, df_charges, all_lines)

Compute bond-projected E-fields using point charges only.

Parameters:

  • bond_indices (list of tuples) – List of (atomA_idx, atomB_idx) tuples specifying bonds.

  • df_charges (pd.DataFrame) – DataFrame with columns: [‘Atom’, ‘x’, ‘y’, ‘z’, ‘charge’].

  • all_lines (list) – Atom indices to include in the E-field calculation.

Returns:

[E_projected, bonded_atoms, bond_indices, bond_lens, E_proj_atomwise, E_proj_atomwise_list]

Return type:

list

vectorEfield(vector_points, xyz_filepath, atom_multipole_file, all_lines, bool_multipole, df_ptchg=None)

Compute E-field projected along vectors defined by pairs of Cartesian points.

Like bondEfield(), but the projection axis is defined by explicit Cartesian coordinates rather than atom indices. E-field is evaluated at both endpoints, averaged, then projected along the A→B unit vector.

Parameters:

  • vector_points (list of tuples) – List of ((x1,y1,z1), (x2,y2,z2)) pairs in Angstrom. Each pair defines a vector from point A to point B.

  • xyz_filepath (str) – Path to XYZ structure file.

  • atom_multipole_file (str) – Path to charge or multipole file.

  • all_lines (list) – Atom indices to include in the calculation.

  • bool_multipole (bool) – If True, use multipole expansion; if False, use monopole charges.

  • df_ptchg (pd.DataFrame or None) – Optional pre-loaded point charge dataframe.

Returns:

[E_projected, vector_labels, vector_points, vector_lens,

E_proj_atomwise, E_proj_atomwise_list]

Return type:

list

vectorEfield_pointcharges(vector_points, df_charges, all_lines)

Compute E-field projections along Cartesian point vectors using point charges only.

Parameters:

  • vector_points (list of tuples) – List of ((x1,y1,z1), (x2,y2,z2)) pairs in Angstrom.

  • df_charges (pd.DataFrame) – DataFrame with columns: [‘Atom’, ‘x’, ‘y’, ‘z’, ‘charge’].

  • all_lines (list) – Atom indices to include in the calculation.

Returns:

[E_projected, vector_labels, vector_points, vector_lens,

E_proj_atomwise, E_proj_atomwise_list]

Return type:

list

esp_coord_shell(metal_idx, charge_file, path_to_xyz, n_shells=1)

Calculate ESP accounting for electrostatic contributions from n coordination shells.

This unified method replaces esp_first_coord() and esp_second_coord() by accepting a parameter to specify how many coordination shells to include.

Parameters:

  • metal_idx (int) – Atom index of the central metal

  • charge_file (str) – Path to charge filename

  • path_to_xyz (str) – Path to xyz structure file

  • n_shells (int, optional) – Number of coordination shells to include (default=1) n_shells=1: first coordination shell only n_shells=2: first and second coordination shells

Returns:

ESP value for the specified coordination shell(s)

Return type:

float

charge_atom(atom_idx)

Compute charges and return total charge and partial charge of atom Input: filename: string of charge filename atom_idx: integer of atom index Output: list of total charge and partial charge of atom

charge_atoms(chg_filename, xyzfilename)

Input: filename: string of charge filename or multipole filename… if the xyzfilename is included must be multipole!

getAtomInfo(atom_idx)

Input: filename: string of charge filename atom_idx: integer of atom index Output: list of total charge and partial charge of atom

getAtomsInfo(atom_indices)

Input: filename: string of charge filename atom_indices: list of integers of atom indices Output: list of total charge and partial charge of atom

esp_bydistance(path_to_xyz, espatom_idx, charge_file)

Input: path_to_xyz: string of xyz filename espatom_idx: integer of atom index charge_file: string of charge filename Output: list of ESP and atomic symbol

getESP(charge_types='CM5', ESPdata_filename='ESP', multiwfn_path='None', use_multipole=False, include_decay=False, include_coord_shells=False, visualize=None, dielectric=1, dielectric_scale=1, includePtChgs=None, ptchg_paths=None, pdb_charge_paths=None)

Unified ESP calculation method supporting multiple modes.

This consolidated method replaces getESPData(), getESPMultipole(), and getESPDecay() to eliminate code duplication while preserving all functionality.

Parameters:

  • charge_types (list of str or str) – Charge partitioning scheme(s). Options: ‘Hirshfeld’, ‘Voronoi’, ‘Mulliken’, ‘Lowdin’, ‘SCPA’, ‘Becke’, ‘ADCH’, ‘CHELPG’, ‘MK’, ‘AIM’, ‘Hirshfeld_I’, ‘CM5’, ‘EEM’, ‘RESP’, ‘PEOE’. Can be a single string (for multipole mode) or list of strings.

  • ESPdata_filename (str) – Output CSV filename (without extension).

  • multiwfn_path (str) – Path to Multiwfn executable. Ignored when pdb_charge_paths is provided.

  • use_multipole (bool, optional) – If True, use multipole expansion; if False, use monopole charges (default: False).

  • include_decay (bool, optional) – If True, include distance-sorted ESP decay analysis (default: False).

  • include_coord_shells (bool, optional) – If True, include first and second coordination shell ESP (default: False).

  • visualize (bool or None, optional) – If True, create PDB file with ESP values in B-factor column. None uses config defaults (default: None).

  • dielectric (float, optional) – Dielectric constant (default: 1).

  • dielectric_scale (float, optional) – Scaling factor for point charges. Point charge contributions are scaled by 1/sqrt(dielectric_scale). Typically used in QM/MM calculations (default: 1).

  • includePtChgs (bool or None, optional) – Override object-level includePtChgs setting for this calculation. If None, uses object-level setting (default: None).

  • ptchg_paths (list of str or None, optional) – Override point charge file paths for this calculation. List of paths (one per job), with None for jobs without point charges. If provided, takes precedence over object-level settings (default: None).

  • pdb_charge_paths (str or list of str or None, optional) – If provided, use partial charges from PDB B-factor columns only (no Multiwfn). Can be a single PDB path or a list of paths (one per job), with None to skip a job. Overrides the object-level pdb_charge_paths if the object was created in pdb_only mode. This mode supports monopole charges only (use_multipole must be False).

Returns:

  • pd.DataFrame – DataFrame with ESP results saved to CSV.

  • Output Files

  • ————

  • **CSV file (*** Auto-generated with format:*)

  • - Monopole (esp_{charge_type}.csv)

  • - Multipole (esp_multipole_{charge_type}.csv)

  • - Examples (esp_Hirshfeld_I.csv, esp_multipole_Hirshfeld_I.csv)

  • - Contains ESP values at metal centers for all jobs

  • - Columns (job name, charge type, ESP values, coordinates, etc.)

  • **PDB file (if visualize=True) (** esp_{charge_type}_{structure}_atom{idx}.pdb)

  • - One file per job and metal atom

  • - B-factor column contains ESP values at each atom

  • - Example (esp_Hirshfeld_I_complex_atom25.pdb)

  • Note (If you provide a custom filename (not ‘output’, ‘esp’, etc.),)

  • your filename will be used instead of automatic naming.

Examples

>>> # Monopole mode (replaces getESPData)
>>> es.getESP(['Hirshfeld', 'Hirshfeld_I'], 'esp_mono', ...)

>>> # Multipole mode (replaces getESPMultipole)
>>> es.getESP('Hirshfeld_I', 'esp_multi', ..., use_multipole=True)

>>> # Decay analysis mode (replaces getESPDecay)
>>> es.getESP(['Hirshfeld'], 'esp_decay', ..., include_decay=True, include_coord_shells=True)

>>> # pdb_only mode (no Multiwfn required)
>>> es.getESP(pdb_charge_paths=['job1.pdb', 'job2.pdb'])
getEfield(charge_types='Hirshfeld_I', Efielddata_filename='ef', multiwfn_path=None, multipole_bool=False, input_bond_indices=[], auto_find_bonds=False, save_atomwise_decomposition=False, visualize=None, dielectric=1, dielectric_scale=1, includePtChgs=None, ptchg_paths=None, pdb_charge_paths=None)

Unified E-field calculation method supporting multiple modes.

This consolidated method replaces getEFieldMultipole(), getEfield_acrossBond(), and getEfield_decomposable() to eliminate code duplication while preserving all functionality.

Parameters:

  • charge_types (str or list of str) – Charge partitioning scheme. If a list is provided, only the first entry is used. Options: ‘Hirshfeld’, ‘Becke’, ‘Hirshfeld_I’, etc. When pdb_charge_paths is provided, this value is used for labeling only.

  • Efielddata_filename (str) – Output CSV filename (without extension).

  • multiwfn_path (str) – Path to Multiwfn executable.

  • multipole_bool (bool, optional) – If True, use multipole expansion; if False, use monopole charges (default: False).

  • input_bond_indices (list, optional) – Bond indices as list of tuples [(atomA, atomB), …]. If empty, bonds are auto-found using metal indices (esp_atom_idx). If provided, bonds are used as-is unless auto_find_bonds=True.

  • auto_find_bonds (bool, optional) – If True, automatically find bonds to adjacent atoms using getBondedAtoms() and ignore input_bond_indices (default: False).

  • save_atomwise_decomposition (bool, optional) – If True, save atom-wise E-field decomposition to CSV (default: False).

  • visualize (bool or None, optional) – Override config visualization settings. None uses config defaults (default: None).

  • dielectric (float, optional) – Dielectric constant (default: 1).

  • dielectric_scale (float, optional) – Scaling factor for point charges. Point charge contributions are scaled by 1/sqrt(dielectric_scale). Typically used in QM/MM calculations (default: 1).

  • includePtChgs (bool or None, optional) – Override object-level includePtChgs setting for this calculation. If None, uses object-level setting (default: None).

  • ptchg_paths (list of str or None, optional) – Override point charge file paths for this calculation. List of paths (one per job), with None for jobs without point charges. If provided, takes precedence over object-level settings (default: None).

  • pdb_charge_paths (str or list of str or None, optional) – If provided, use partial charges from PDB B-factor columns only (no Multiwfn). Can be a single PDB path or a list of paths (one per job), with None to skip a job. This mode supports monopole charges only (multipole_bool must be False).

Returns:

  • pd.DataFrame or tuple of pd.DataFrame – If save_atomwise_decomposition=False: returns DataFrame with E-field results If save_atomwise_decomposition=True: returns (total_df, atomwise_df) tuple

  • Output Files

  • ————

  • **Main CSV file (*** Auto-generated with format:*)

  • - Monopole (ef_{charge_type}.csv)

  • - Multipole (ef_multipole_{charge_type}.csv)

  • - Examples (ef_Hirshfeld_I.csv, ef_multipole_Hirshfeld_I.csv)

  • - Contains E-field magnitudes and projections for all jobs and bonds

  • - Columns (job name, bond atoms, E-field magnitude (MV/cm), projections, etc.)

  • **Atomwise CSV (if save_atomwise_decomposition=True) (****)

  • - Format (ef_{charge_type}_atomwise.csv or ef_multipole_{charge_type}_atomwise.csv)

  • - Example (ef_Hirshfeld_I_atomwise.csv)

  • - Per-atom E-field contributions for each bond

  • - Shows which atoms contribute most to the E-field

  • - Columns (job, bond, atom index, contribution magnitude, etc.)

  • **PDB files (if visualize=True) (****)

  • - Per bond (ef_{charge_type}_{structure}_bond{i}-{j}.pdb)

  • - Combined (ef_{charge_type}_{structure}_combined.pdb)

  • - B-factor contains E-field contribution from each atom

  • - Examples (ef_Hirshfeld_I_complex_bond25-26.pdb)

  • Note (If you provide a custom filename (not ‘output’, ‘efield’, etc.),)

  • your filename will be used instead of automatic naming.

Examples

>>> # Basic multipole mode with specified bonds
>>> es.getEfield('Hirshfeld_I', 'efield', ..., input_bond_indices=[(0,1), (0,2)])

>>> # Auto-find bonds to adjacent atoms
>>> es.getEfield('Hirshfeld_I', 'efield', ..., auto_find_bonds=True)

>>> # Monopole mode with atom-wise decomposition
>>> es.getEfield('Hirshfeld', 'efield', ..., multipole_bool=False, save_atomwise_decomposition=True)
getEfield_byResidue(residue_groups, input_bond_indices, exclude_atoms=None, dielectric=1, pdb_charge_paths=None, output_filename='ef_byresidue', other_label='Other')

Decompose bond-projected E-field contributions by PDB residue name.

Computes the same bond-projected E-field as getEfield(…, pdb_charge_paths=…), but instead of returning one value per bond, sums the atom-wise contributions within each user-defined group of residues (e.g. metal centers, ligands, counterions, solvent) so the source of the field can be attributed.

This method only supports PDB B-factor charges (pdb_only mode), since residue names are read from the PDB ‘ResName’ field.

Parameters:

  • residue_groups (dict) – {group_name: [resname, resname, …]}. Maps each group label to the list of PDB residue names (as written in the PDB ResName field) that belong to it. Example: {‘Metal’: [‘PD1’, ‘PD2’], ‘Solvent’: [‘MOL’]}.

  • input_bond_indices (list of list of tuple) – Bond indices per job: input_bond_indices[j] is a list of (atomA_idx, atomB_idx) tuples (0-indexed) for job j.

  • exclude_atoms (list of int, or list of list of int, optional) – Atom indices to exclude entirely from the E-field calculation (e.g. the solute/guest whose own field contribution shouldn’t count). A flat list applies to every job; a list of lists gives per-job exclusions (default: None).

  • dielectric (float, optional) – Dielectric constant (default: 1).

  • pdb_charge_paths (str or list of str or None, optional) – PDB file path(s) with charges in the B-factor column, one per job. If None, uses the object-level pdb_charge_paths (pdb_only mode).

  • output_filename (str, optional) – Output CSV filename without extension (default: ‘ef_byresidue’).

  • other_label (str or None, optional) – Label for atoms whose residue name isn’t listed in residue_groups. If None, unmatched atoms are dropped from the output (default: ‘Other’).

Returns:

One row per (Structure, Bond, Group) with columns: [‘Structure’, ‘Bond’, ‘AtomA_idx’, ‘AtomB_idx’, ‘Group’, ‘N_atoms’, ‘Efield_V_per_A’, ‘Efield_MV_per_cm’].

Return type:

pd.DataFrame

Examples

>>> es = Electrostatics(pdb_charge_paths=['frame1.pdb', 'frame2.pdb'])
>>> df = es.getEfield_byResidue(
...     residue_groups={'Metal': ['PD1', 'PD2'], 'Solvent': ['MOL']},
...     input_bond_indices=[[(10, 11)], [(10, 11)]],
...     exclude_atoms=[0, 1, 2],
... )
getESP_byResidue(residue_groups, esp_atom_idx=None, exclude_atoms=None, dielectric=1, pdb_charge_paths=None, output_filename='esp_byresidue', other_label='Other')

Decompose ESP contributions at an atom by PDB residue name.

Computes the same monopole ESP as getESP(…, pdb_charge_paths=…) at each specified atom, but sums the atom-wise contributions within each user-defined group of residues (e.g. metal centers, ligands, counterions, solvent) so the source of the potential can be attributed.

This method only supports PDB B-factor charges (pdb_only mode), since residue names are read from the PDB ‘ResName’ field.

Parameters:

  • residue_groups (dict) – {group_name: [resname, resname, …]}. Maps each group label to the list of PDB residue names (as written in the PDB ResName field) that belong to it. Example: {‘Metal’: [‘PD1’, ‘PD2’], ‘Solvent’: [‘MOL’]}.

  • esp_atom_idx (list of int, or None, optional) – Atom index (0-indexed) at which to evaluate the ESP, one per job. If None, uses the object-level esp_atom_idx provided at construction (default: None).

  • exclude_atoms (list of int, or list of list of int, optional) – Atom indices to exclude entirely from the ESP calculation. A flat list applies to every job; a list of lists gives per-job exclusions (default: None).

  • dielectric (float, optional) – Dielectric constant (default: 1).

  • pdb_charge_paths (str or list of str or None, optional) – PDB file path(s) with charges in the B-factor column, one per job. If None, uses the object-level pdb_charge_paths (pdb_only mode).

  • output_filename (str, optional) – Output CSV filename without extension (default: ‘esp_byresidue’).

  • other_label (str or None, optional) – Label for atoms whose residue name isn’t listed in residue_groups. If None, unmatched atoms are dropped from the output (default: ‘Other’).

Returns:

One row per (Structure, Group) with columns: [‘Structure’, ‘ESP_Atom_idx’, ‘ESP_Atom_symbol’, ‘Group’, ‘N_atoms’, ‘ESP_V’].

Return type:

pd.DataFrame

Examples

>>> es = Electrostatics(pdb_charge_paths=['frame1.pdb', 'frame2.pdb'])
>>> df = es.getESP_byResidue(
...     residue_groups={'Metal': ['PD1', 'PD2'], 'Solvent': ['MOL']},
...     esp_atom_idx=[12, 12],
...     exclude_atoms=[0, 1, 2],
... )
getEfield_vectorProjection(charge_types='Hirshfeld_I', Efielddata_filename='ef_vec', multiwfn_path=None, multipole_bool=False, input_vector_points=[], save_atomwise_decomposition=False, dielectric=1, dielectric_scale=1, includePtChgs=None, ptchg_paths=None, pdb_charge_paths=None)

Compute E-field projected along user-defined vectors in Cartesian space.

Like getEfield(), but the projection axis is defined by explicit Cartesian coordinates instead of atom pairs. This allows computing E-field projections along any arbitrary vector in space, not just along chemical bonds.

The E-field is evaluated at both endpoints of each vector, averaged, then projected along the A→B unit vector — matching the convention of getEfield().

Parameters:

  • charge_types (str or list of str) – Charge partitioning scheme. If a list is provided, only the first is used.

  • Efielddata_filename (str) – Output CSV filename (without extension).

  • multiwfn_path (str) – Path to Multiwfn executable.

  • multipole_bool (bool, optional) – If True, use multipole expansion; if False, use monopole charges (default: False).

  • input_vector_points (list) –

    One entry per job, each entry is a list of ((x1,y1,z1), (x2,y2,z2)) tuples (in Angstrom) defining the projection vectors for that job. Example for two jobs with one vector each:

    [
    [((0.0, 0.0, 0.0), (1.0, 0.0, 0.0))],   # job 0
    [((0.5, 0.5, 0.0), (0.5, 0.5, 1.0))],   # job 1
]

  • save_atomwise_decomposition (bool, optional) – If True, save per-atom E-field contributions to CSV (default: False).

  • dielectric (float, optional) – Dielectric constant (default: 1).

  • dielectric_scale (float, optional) – Scaling factor for point charges (default: 1).

  • includePtChgs (bool or None, optional) – Override object-level includePtChgs setting. None uses object-level setting.

  • ptchg_paths (list of str or None, optional) – Override point charge file paths (one per job).

  • pdb_charge_paths (str or list of str or None, optional) – If provided, use partial charges from PDB B-factor columns only (no Multiwfn).

Returns:

If save_atomwise_decomposition=False: DataFrame with E-field results. If save_atomwise_decomposition=True: (total_df, atomwise_df) tuple.

Return type:

pd.DataFrame or tuple of pd.DataFrame

Examples

>>> # Project along the x-axis at the origin for all jobs
>>> vectors = [((0.0, 0.0, 0.0), (1.0, 0.0, 0.0))]
>>> ef_df = estat.getEfield_vectorProjection(
...     'Hirshfeld_I', 'ef_vec', multiwfn_path,
...     input_vector_points=[vectors, vectors]
... )
getDipoleDifference(charge_types='Hirshfeld_I', multiwfn_path=None, output_filename='dipole_diff', align_heavy_only=True, align_indices=None, pdb_charge_paths=None, save_aligned_xyz=True)

Compute the difference dipole vector between two molecular structures.

Kabsch-aligns the second structure onto the first (minimizing RMSD), then computes the partial-charge dipole vector for both in the common aligned frame. The difference dipole Δμ = μ₂ − μ₁ describes the direction in which the charge distribution shifts between the two states.

Requires exactly two jobs in the Electrostatics object (one per structure).

Parameters:

  • charge_types (str or list of str, optional) – Charge partitioning scheme. If a list, only the first is used. Default: ‘Hirshfeld_I’.

  • multiwfn_path (str or None) – Path to Multiwfn executable. Not needed when pdb_charge_paths is set.

  • output_filename (str, optional) – Base name for output CSV (default: ‘dipole_diff’).

  • align_heavy_only (bool, optional) – If True (default), use only non-hydrogen atoms for Kabsch alignment. Ignored when align_indices is specified.

  • align_indices (array-like or None, optional) – Explicit 0-based atom indices to use for alignment. Overrides align_heavy_only when provided.

  • pdb_charge_paths (list of str or None, optional) – If provided, read partial charges from PDB B-factor columns instead of running Multiwfn. Must be a list of two paths, one per job.

  • save_aligned_xyz (bool, optional) – If True (default), save the aligned structure 2 as an XYZ file in structure 1’s coordinate frame (useful for visualization).

Returns:

  • dict with keys – ‘dipole_1’ : np.ndarray (3,), μ₁ in Debye ‘dipole_2_aligned’ : np.ndarray (3,), μ₂ in Debye (aligned to struct 1) ‘delta_dipole’ : np.ndarray (3,), Δμ = μ₂ − μ₁ in Debye ‘delta_dipole_magnitude’ : float, |Δμ| in Debye ‘delta_dipole_direction’ : np.ndarray (3,), unit vector of Δμ ‘rotation_matrix’ : np.ndarray (3,3), Kabsch rotation (struct 2 → struct 1) ‘rmsd’ : float, RMSD after alignment in Angstrom ‘charge_type’ : str ‘xyz_1’ : str ‘xyz_2’ : str

  • Also saves –

    • {output_filename}.csv — scalar results

    • {output_filename}_struct2_aligned.xyz (if save_aligned_xyz=True)

Notes

Dipoles are computed relative to each structure’s geometric centroid, making them origin-independent for neutral molecules. For charged systems the dipole is origin-dependent; the centroid convention is used throughout.

Examples

>>> estat = Electrostatics(
...     molden_paths=['state1.molden', 'state2.molden'],
...     xyz_paths=['state1.xyz', 'state2.xyz']
... )
>>> result = estat.getDipoleDifference('Hirshfeld_I', '/path/to/multiwfn')
>>> print(result['delta_dipole_direction'])   # unit vector of charge shift
>>> print(result['delta_dipole_magnitude'])   # |Δμ| in Debye
getCharges(charge_types='Hirshfeld_I', multiwfn_path=None, output_filename='charges', write_pdb=False, pdb_bfactor=True)

Compute partial charges for all atoms using specified charge partitioning scheme(s).

This method provides a standalone interface for computing partial charges without calculating E-fields or ESP. Charges can optionally be written to PDB files with the B-factor column containing the charge values.

Parameters:

  • charge_types (str or list of str) – Charge partitioning scheme(s) to use. Options include: ‘RESP’, ‘Hirshfeld’, ‘Hirshfeld_I’, ‘CHELPG’, ‘Mulliken’, ‘Lowdin’, ‘SCPA’, ‘Becke’, ‘ADCH’, ‘Voronoi’, ‘MK’, ‘AIM’, ‘CM5’, ‘EEM’, ‘PEOE’

  • multiwfn_path (str) – Path to the Multiwfn executable.

  • output_filename (str, optional) – Base name for output files (default: ‘charges’). CSV will be saved as ‘{output_filename}.csv’.

  • write_pdb (bool, optional) – If True, create PDB files for each structure (default: False).

  • pdb_bfactor (bool, optional) – If True and write_pdb=True, store charges in B-factor column (default: True). This allows visualization of charges in molecular viewers like PyMOL.

Returns:

  • pd.DataFrame – DataFrame containing partial charges for all atoms in each structure. Columns include: ‘Job’, ‘Charge_Type’, ‘Atom_Index’, ‘Element’, ‘x’, ‘y’, ‘z’, ‘Charge’.

  • Output Files

  • ————

  • **CSV file (** {output_filename}.csv) – Contains all partial charges for all atoms and all jobs.

  • **PDB files (if write_pdb=True) (** {structure}_{charge_type}_charges.pdb) – One PDB file per job per charge type. B-factor column contains the partial charge values.

Examples

>>> es = Electrostatics(molden_paths, xyz_paths)
>>>
>>> # Get RESP charges for all atoms
>>> df = es.getCharges('RESP', '/path/to/multiwfn')
>>>
>>> # Get multiple charge types with PDB output
>>> df = es.getCharges(['RESP', 'Hirshfeld_I'], '/path/to/multiwfn',
...                    output_filename='my_charges', write_pdb=True)
>>>
>>> # Access charges for a specific atom
>>> atom_5_charges = df[df['Atom_Index'] == 5]

Notes

  • Charges are computed via Multiwfn and cached in ‘Charges{type}.txt’ files.

  • Subsequent calls with the same charge type will reuse cached results unless config[‘rerun’] is True.

  • PDB files with B-factor charges can be visualized in PyMOL using: spectrum b, blue_white_red, minimum=-1, maximum=1

get_residues(auto_gen, solute_solvent_dict)
partitionCharge(multipole_bool, f, multiwfn_path, charge_type, owd)

Partition electron density using Multiwfn to generate partial charges or multipole moments.

This is the centralized interface for all charge partitioning schemes (Hirshfeld, CHELPG, etc.). Checks if calculation was previously completed and reuses results unless rerun=True.

Parameters:

multipole_boolbool

True to compute multipole moments, False for monopoles only

fstr

Complete path to folder containing the calculation

multiwfn_pathstr

Path to Multiwfn executable

charge_typestr

Partitioning scheme (e.g., ‘Hirshfeld’, ‘CHELPG’, ‘Hirshfeld_I’)

owdstr

Original working directory

Returns:

comp_costfloat

Computation time in seconds, 0 if calculation was previously completed (skipped), or -1 if calculation failed

Notes:

For Hirshfeld-I calculations, atmrad files are automatically loaded from pyef.resources.atmrad (bundled with the package).

get_residueDipoles(charge_type, multiwfn_path, multipole_bool, solute_indices=[], num_atoms_solvent=0, save_distance_dipoles=False)

Function computes a partial charges of all atoms in one system… if the desired file already exists the just return it!! Inputs: ——- charge_type: string, possible choices include: ‘Hirshfeld’, ‘Voronoi’, ‘Mulliken’, ‘Lowdin’, ‘SCPA’, ‘Becke’, ‘ADCH’, ‘CHELPG’, ‘MK’, ‘AIM’, ‘Hirshfeld_I’, ‘CM5’, ‘EEM’, ‘RESP’, ‘PEOE’}

Also accepts multipole options including: ‘Hirshfeld’: [‘3’, ‘2’], ‘Hirshfeld_I’: [‘4’, ‘1’, ‘2’], ‘Becke’: [‘1’, ‘2’]

getpartialchgs(charge_types, lst_atom_idxs, partial_chg_filename, multiwfn_path)

Function computes partial charges on a select set of atoms using the charge scheme specified in charge types. Note atom indices will be carried over between csvs

Inputs:

– charge_types: list(str)

list of strings

lst_of_atom_idxs: list(int)

list of integers denoting atom indices (0 indexed!)

partial_chg_filename: string

Name of the output file name

multiwfn_path: string

Path to the multiwfn executable

Path to the atmrad executable

Outputs:

– df: Pandas Dataframe with partial charge info

Will Create a csv file entitled partial_chg_filename.csv with partial charge info

compute_dipole(res_chgs, nuc_chgs)

res_coords in units of angstroms, res_chgs are in a.u. units returns dipole in units of Debye

getdipole_residues(res_dict, df)

res_dict: dictionary with strings mapped to list of lists(int) strings denote name of residues which are mapped to the associated atom indices in the lists(0 indexed!)

getcharge_residues(charge_types, res_dict, partial_chg_filename, multiwfn_path, multipole_mode=True, polarization_scheme='Hirshfeld_I')

Function computes partial charges on a select set of atoms using the charge scheme specified in charge types. Note atom indices will be carried over between csvs

Inputs:

charge_types: list(str)

list of strings

res_dict: dictionary with strings mapped to list of lists(int)

strings denote name of residues which are mapped to the associated atom indices in the lists(0 indexed!)

partial_chg_filename: string

Name of the output file name

multiwfn_path: string

Path to the multiwfn executable

Path to the atmrad executable

multipole_mode: boolean

If True, will use multipole mode to compute partial charges, if False will use the standard method

Outputs:

Notes

Will Create a csv file entitled partial_chg_filename.csv with partial charge info

compute_atomwise_dipole_contributions(atomwise_Efield, df_substrate, df_env, df_all_atoms, one_mol, filename)

Compute atom-wise energy contributions from environment E-fields interacting with substrate dipoles.

The interaction energy is: E = -μ·E (negative dot product of dipole moment with electric field)

Parameters:

atomwise_Efielddict

Dictionary with (substrate_idx, env_idx) as keys and E-field vector [Ex, Ey, Ez] (in V/Å) as values

df_substrateDataFrame

DataFrame containing substrate atom information

df_envDataFrame

DataFrame containing environment atom information

df_all_atomsDataFrame

DataFrame containing all atom information (for element names)

one_molfloat

Avogadro’s number

filenamestr

Name of the current file being processed

Returns:

DataFrame with columns: FileName, Env_Atom_Index, Env_Element, Total_Energy_Contribution_kcal_mol,

Substrate_Atom_Index, Substrate_Element, Energy_Contribution_kcal_mol

visualize_atomwise_contributions(df_atomwise, file_path_xyz, charge_file, filename, charge_type)

Create PDB visualization of atom-wise energy contributions.

Parameters:

df_atomwiseDataFrame

DataFrame with atom-wise contribution data from _compute_atomwise_contributions

file_path_xyzstr

Path to XYZ file

charge_filestr

Path to charge/multipole file

filenamestr

Folder/file name for output

charge_typestr

Type of charge calculation used

compute_atomwise_contributions(atomwise_ESP, df_substrate, df_env, df_all_atoms, one_mol, filename)

Compute atom-wise energy contributions from environment atoms to substrate stabilization.

Parameters:

atomwise_ESPdict

Dictionary with (substrate_idx, env_idx) as keys and ESP contribution (in Volts) as values

df_substrateDataFrame

DataFrame containing substrate atom information

df_envDataFrame

DataFrame containing environment atom information

df_all_atomsDataFrame

DataFrame containing all atom information (for element names)

one_molfloat

Avogadro’s number

filenamestr

Name of the current file being processed

Returns:

DataFrame with columns: FileName, Env_Atom_Index, Env_Element, Total_Energy_Contribution_kcal_mol,

Substrate_Atom_Index, Substrate_Element, Energy_Contribution_kcal_mol

build_interaction_tensor(r_i, r_j, multipole_order, dielectric)

Build the multipole interaction tensor T_ij using AMOEBA formalism.

This tensor contains derivatives of 1/r and enables computation of all multipole-multipole interactions via: E = M_i^T · T_ij · M_j

Parameters:

r_i, r_jarray-like

Position vectors of atoms i and j in Angstroms

multipole_orderint

Maximum multipole order (1=monopole, 2=+dipole, 3=+quadrupole)

dielectricfloat

Dielectric constant

Returns:

Tndarray

Interaction tensor with dimensions based on multipole_order - Order 1: 1×1 (monopole-monopole only) - Order 2: 4×4 (monopole + dipole) - Order 3: 9×9 (monopole + dipole + quadrupole, traceless)

Notes:

Tensor elements represent spatial derivatives of Coulomb potential: - T[0,0]: 1/r (monopole-monopole) - T[0,1:4]: ∂(1/r)/∂r_j (monopole-dipole) - T[1:4,1:4]: ∂²(1/r)/∂r_i∂r_j (dipole-dipole) - Higher orders follow similar pattern

build_multipole_vector(atom_data, multipole_order)

Build multipole moment vector M for an atom using AMOEBA formalism.

Parameters:

atom_datadict or pd.Series

Atom data containing ‘Atom_Charge’, ‘Dipole_Moment’, ‘Quadrupole_Moment’

multipole_orderint

Maximum multipole order (1=monopole, 2=+dipole, 3=+quadrupole)

Returns:

Mndarray

Multipole vector [q, μ_x, μ_y, μ_z, Q_xx, Q_xy, …] in SI-based units - Order 1: [q] (1 element) - Order 2: [q, μ_x, μ_y, μ_z] (4 elements) - Order 3: [q, μ_x, μ_y, μ_z, Q_xx, Q_xy, Q_xz, Q_yy, Q_yz] (9 elements, traceless)

Notes:

  • Charges are dimensionless (in units of elementary charge)

  • Dipoles are converted from Bohr to meters

  • Quadrupoles are converted from Bohr² to meters²

  • Uses traceless quadrupole representation (Q_zz = -Q_xx - Q_yy)

getElectrostatic_stabilization(multiwfn_path=None, substrate_idxs=[], charge_type='Hirshfeld_I', name_dataStorage='estatic', env_idxs=None, save_atomwise_decomposition=False, visualize=None, multipole_order=2, substrate_multipole_order=None, env_multipole_order=None, dielectric=1, dielectric_scale=1, includePtChgs=None, ptchg_paths=None, pdb_charge_paths=None)

Compute electrostatic stabilization using AMOEBA-style multipole tensor formalism.

This function implements the complete multipole expansion for intermolecular electrostatic interactions using the tensor formalism:

E_total = Σᵢ∈substrate Σⱼ∈environment M_i^T · T_ij · M_j

where M is the multipole vector and T_ij is the interaction tensor containing all derivatives of the Coulomb potential. This naturally includes ALL multipole-multipole interaction terms:

  • Order 1 (monopole): charge-charge (q×q)

  • Order 2 (+ dipole): + charge-dipole (q×μ) + dipole-dipole (μ×μ)

  • Order 3 (+ quadrupole): + charge-quadrupole (q×Q) + dipole-quadrupole (μ×Q) + quadrupole-quadrupole (Q×Q)

This is the formalism used in the AMOEBA polarizable force field and provides a complete treatment of distributed multipole interactions.

Parameters:

multiwfn_pathstr

Path to Multiwfn executable

substrate_idxslist

List of substrate atom indices (the molecule being stabilized)

charge_typestr, optional

Charge partitioning scheme (default: ‘Hirshfeld_I’)

name_dataStoragestr, optional

Output filename prefix (default: ‘estaticFile_tensor’)

env_idxslist or None, optional

List of environment atom indices. If None, uses all non-substrate atoms (default: None)

save_atomwise_decompositionbool, optional

If True, save atom-wise decomposition to CSV showing each substrate atom’s contribution (default: False)

visualizebool or None, optional

If True, create PDB files with energy contributions in B-factor column. None uses config defaults (default: None)

multipole_orderint, optional

Order of multipole expansion: - 1: monopole only (charge-charge) - 2: monopole + dipole (includes charge-charge, charge-dipole, dipole-dipole) - 3: monopole + dipole + quadrupole (includes all terms up to Q×Q) (default: 2)

substrate_multipole_orderint or None, optional

Multipole order for substrate atoms. If None, uses multipole_order. Set to 1 for QM/MM where substrate is MM (charges only). (default: None)

env_multipole_orderint or None, optional

Multipole order for environment atoms. If None, uses multipole_order. Set to 1 for QM/MM where environment is MM (charges only). (default: None)

dielectricfloat, optional

Dielectric constant (default: 1).

dielectric_scalefloat, optional

Scaling factor for point charges. Point charge contributions are scaled by 1/sqrt(dielectric_scale). Typically used in QM/MM calculations (default: 1).

includePtChgsbool or None, optional

Include point charges as environment atoms only. If None, uses object-level setting (default: None). Point charges are ONLY added to environment, never to substrate.

ptchg_pathslist of str or None, optional

Override point charge file paths for this calculation. List of paths (one per job), with None for jobs without point charges. Point charges are treated as monopole-only environment atoms. If provided, takes precedence over object-level settings (default: None).

pdb_charge_pathsstr or list of str or None, optional

If provided, use partial charges from PDB B-factor columns only (no Multiwfn). Monopole-only (multipole_order is forced to 1). Coordinates are also read from the PDB, so no separate xyz file is needed. Overrides the object-level pdb_charge_paths when the object was created in pdb_only mode.

Returns:

pd.DataFrame or tuple of pd.DataFrame

If save_atomwise_decomposition=False: returns DataFrame with total stabilization energies If save_atomwise_decomposition=True: returns (total_df, atomwise_df) tuple

Output Files

Main CSV file: Auto-generated with format: - Monopole: estab_{charge_type}.csv - Multipole: estab_multipole{N}_{charge_type}.csv (N = multipole order) - Examples: estab_Hirshfeld_I.csv, estab_multipole2_Hirshfeld_I.csv - Contains total electrostatic stabilization energies (kcal/mol) for all jobs - Columns: job name, total energy, substrate/environment details, etc.

Atomwise CSV (if save_atomwise_decomposition=True): - Format: estab_{charge_type}_atomwise.csv or estab_multipole{N}_{charge_type}_atomwise.csv - Examples: estab_Hirshfeld_I_atomwise.csv, estab_multipole2_Hirshfeld_I_atomwise.csv - Per-atom energy contributions showing which atoms contribute most - Columns: job, atom index, energy contribution, coordinates, etc.

PDB file (if visualize=True): estab_{charge_type}_{structure}_tensor_sub{N}_env{M}.pdb - B-factor contains energy contribution from each environment atom to substrate - N = substrate multipole order, M = environment multipole order - Example: estab_Hirshfeld_I_complex_tensor_sub2_env1.pdb

Note: If you provide a custom filename (not ‘output’, ‘estatic’, etc.), your filename will be used instead of automatic naming.

Notes:

Comparison with other methods: - getElectrostatic_stabilization(): Computes q·V (substrate charge in environment potential) - get_Electrostatic_stabilization_dipole(): Computes -μ·E (substrate dipole in environment field) - This function: Computes DIRECT pairwise multipole-multipole interactions

Advantages of tensor formalism: - Complete: automatically includes all multipole-multipole terms - Symmetric: properly treats both substrate and environment multipoles - Efficient: single matrix multiplication per atom pair - Extensible: easy to add higher-order terms

Important: - Requires Multiwfn multipole analysis for orders ≥ 2 - Automatic fallback: If multipole files don’t exist, automatically uses charges-only files - QM/MM support: Substrate and environment can have different multipole orders (e.g., QM substrate with multipoles, MM environment with charges only) - Results are the electrostatic STABILIZATION energy (environment → substrate) - Energy is in kcal/mol - Follows AMOEBA force field formalism (Ren & Ponder, J. Phys. Chem. B 2003)

Examples:

>>> # Compute monopole + dipole + dipole-dipole interactions
>>> df = estat.get_Electrostatic_stabilization_tensor(
...     multiwfn_path,
...     substrate_idxs=[[0, 1, 2]],  # 3-atom substrate
...     multipole_order=2  # Include dipole-dipole terms
... )

>>> # Include quadrupole terms with atom-wise decomposition
>>> df, df_atomwise = estat.get_Electrostatic_stabilization_tensor(
...     multiwfn_path,
...     substrate_idxs=[[0, 1, 2]],
...     multipole_order=3,
...     save_atomwise_decomposition=True
... )

>>> # QM/MM: QM substrate (with dipoles) interacting with MM environment (charges only)
>>> df = estat.get_Electrostatic_stabilization_tensor(
...     multiwfn_path,
...     substrate_idxs=[[0, 1, 2]],  # QM region
...     substrate_multipole_order=2,  # QM: charges + dipoles
...     env_multipole_order=1         # MM: charges only
... )
getEfieldsNormalMode(fchk_path, mode_index, output_filename, multiwfn_path, n_solute, charge_type='CM5', mode_vectors=None, align_heavy_only=True, dielectric_scale=1.0)

Project the environmental electric field onto a vibrational normal mode.

Computes the Vibrational Stark Effect (VSE) field:

F_mode = Σ_i F_i · (R @ L_i)

where F_i is the Coulomb field at solute atom i from the environment (QM solvent charges + MM point charges if QM/MM), L_i is the normal mode displacement vector, and R is the Kabsch rotation that aligns the QM reference geometry (from fchk) to each frame’s optimised geometry.

Supports three calculation modes depending on how the Electrostatics object was constructed:

  • Pure QM — only molden files, no point charges

  • QM/MM — molden files + MM point charges (via ptchg_paths / includePtChgs)

  • Reuse existing charges — if Charges{charge_type}.txt already exists in each job folder (e.g. from a previous getEfield run), Multiwfn is skipped for that frame.

Parameters:

  • fchk_path (str) –

    Path to the .fchk file.

    • Gaussian format: geometry + normal modes are read directly; mode_vectors is ignored.

    • Q-Chem format: only reference geometry is read; you must supply mode_vectors.

  • mode_index (int) – 1-based mode number (as labelled by Q-Chem / Gaussian output).

  • output_filename (str) – Stem for the output CSV (extension .csv is added if absent).

  • multiwfn_path (str) – Path to the Multiwfn executable.

  • n_solute (int) – Number of solute atoms (first n_solute rows of each charge file are treated as the molecule of interest; the rest are environment).

  • charge_type (str, optional) – Charge partitioning scheme (default 'CM5').

  • mode_vectors (np.ndarray or None, optional) – Normal mode displacement array, shape (n_modes, n_atoms, 3) or (n_atoms, 3) for a single mode. Required for Q-Chem fchk; ignored for Gaussian fchk.

  • align_heavy_only (bool, optional) – If True (default), use only non-hydrogen atoms for the Kabsch alignment (more stable; avoids H-atom positional noise).

  • dielectric_scale (float, optional) – Scale factor applied to MM point charges (default 1.0). Set to 1/epsilon for implicit-solvent screening of the MM region.

Returns:

  • pd.DataFrame – Columns: job, F_mode_VA (V/Å), F_mode_MV_cm (MV/cm = 100 × V/Å), plus freq_cm1 if read from Gaussian fchk.

  • Output Files

  • ————

  • {output_filename}.csv — one row per job.

Examples

Pure QM (Gaussian fchk, mode from fchk):

>>> es = Electrostatics(molden_paths, xyz_paths)
>>> es.prepData()
>>> df = es.getEfieldsNormalMode(
...     fchk_path='freq.fchk',
...     mode_index=33,
...     output_filename='nm_ef_CM5',
...     multiwfn_path='/path/to/Multiwfn',
...     n_solute=15,
... )

QM/MM (Q-Chem fchk, mode vectors supplied):

>>> es = Electrostatics(molden_paths, xyz_paths, ptchg_paths=ptchg_list)
>>> es.prepData()
>>> df = es.getEfieldsNormalMode(
...     fchk_path='freq.fchk',
...     mode_index=33,
...     output_filename='nm_ef_QMMMmode33',
...     multiwfn_path='/path/to/Multiwfn',
...     n_solute=15,
...     mode_vectors=PZQ_MODE33,   # (n_atoms, 3) array
... )
visualize_charges_pdb(xyz_file, charges, pdb_name)

Create a PDB file with charges visualized in the b-factor column.

Parameters:

  • xyz_file (str) – Path to the XYZ structure file

  • charges (array-like) – Array of atomic charges to visualize

  • pdb_name (str) – Output PDB filename

Returns:

Writes PDB file to disk

Return type:

None