Mod:Hunt Research Group/dipole moment atom coordinate

Code to Calculate the Dipole Moment from the Center of Mass (CoM), any Atom, or a given set of coordinates and produce xyz file with the dipole vectors

What the program does?

1. The script needs numpy and contextlib packages to run

2. The script can read a Gausslain log file, extract its coordinates, and calculate the center of mass (CoM) from those coordinates

3. Then it extracts ESP and NBO charges from the log file, determines the center of mass, and calculates dipole w.r.t the center of CoM

4. It saves multiple xyz files to the current directory with a dummy atom (X) with CoM coordinates, and dummy atoms X and Bq to show the dipole vector (from X to Bq)

5. It also prints some information on the screen, including the dipole moment calculated by Gaussian (same info is also exported as a .txt file)

Running the python Script

• Running the script is very simple. Simply save the script as calculate_dipole_wrt_atom_or_coordinates_g16.py and run it as follows (considering your molecule file is test.log)

python3 calculate_dipole_wrt_atom_or_coordinates_g16.py test.log [-atom N] [-origin X Y Z]

• -atom can be used to specify the atom number N
• -origin can be used to specify the X Y Z coordinates as the origin of the dipole vector
• if no argument is provided, it will use CoM by default as the dipole origin
• It will produce the following files with their explanation in parentheses.
  • test_ESP_dipole.xyz (molecule with dummy atom (X) as center of vector tail and Bq to show dipole vector arrow w.r.t ESP)
  • test_NBO_dipole.xyz (molecule with dummy atom (X) as center of vector tail and Bq to show dipole vector arrow w.r.t NBO)

Below is the Script

#!/usr/bin/env python3
"""
Script to calculate the center of mass (CoM), center of charge (CoP for ESP and NBO charges) of a molecule from a log file,
then calculate the dipole moment from different charges, and write a new XYZ file with the CoM added as a pseudo-atom 'X'.
Script by Dr. Muhammad Ali Hashmi with the help of ChatGPT (29/10/2025)
Features:
1. Orientation fallback: Standard → Z-Matrix → Input
2. Option to define dipole origin using specific atom (-atom N, 1-based)
3. Option to define manual origin (-origin x y z)
4. Generates .xyz files with dummy atoms (X, Bq) representing dipole direction
   - X = origin
   - Bq = 6 Å along dipole vector
   - Appended to original atoms
   - Separate files for ESP and NBO dipoles
"""

import numpy as np
import os
import sys
from contextlib import redirect_stdout

# === Atomic weights dictionary ===
atomic_masses = {
"H" : 1.00794, "He": 4.002602, "Li": 6.941, "Be": 9.012182, "B": 10.811,
"C": 12.0107, "N": 14.0067, "O": 15.9994, "F": 18.9984032, "Ne": 20.1797,
"Na": 22.98976928, "Mg": 24.3050, "Al": 26.9815386, "Si": 28.0855,
"P": 30.973762, "S": 32.065, "Cl": 35.453, "Ar": 39.948, "K": 39.0983,
"Ca": 40.078, "Sc": 44.955912, "Ti": 47.867, "V": 50.9415, "Cr": 51.9961,
"Mn": 54.938045, "Fe": 55.845, "Co": 58.933195, "Ni": 58.6934, "Cu": 63.546,
"Zn": 65.38, "Ga": 69.723, "Ge": 72.64, "As": 74.92160, "Se": 78.96,
"Br": 79.904, "Kr": 83.798, "Rb": 85.4678, "Sr": 87.62, "Y": 88.90585,
"Zr": 91.224, "Nb": 92.90638, "Mo": 95.96, "Tc": 98.0, "Ru": 101.07,
"Rh": 102.90550, "Pd": 106.42, "Ag": 107.8682, "Cd": 112.411, "In": 114.818,
"Sn": 118.710, "Sb": 121.760, "Te": 127.60, "I": 126.90447, "Xe": 131.293,
"Cs": 132.9054519, "Ba": 137.327, "La": 138.90547, "Ce": 140.116,
"Pr": 140.90765, "Nd": 144.242, "Pm": 145.0, "Sm": 150.36, "Eu": 151.964,
"Gd": 157.25, "Tb": 158.92535, "Dy": 162.500, "Ho": 164.93032, "Er": 167.259,
"Tm": 168.93421, "Yb": 173.054, "Lu": 174.9668, "Hf": 178.49, "Ta": 180.94788,
"W": 183.84, "Re": 186.207, "Os": 190.23, "Ir": 192.217, "Pt": 195.084,
"Au": 196.966569, "Hg": 200.59, "Tl": 204.3833, "Pb": 207.2, "Bi": 208.98040,
"Po": 209.0, "At": 210.0, "Rn": 222.0, "Fr": 223.0, "Ra": 226.0, "Ac": 227.0,
"Th": 232.03806, "Pa": 231.03588, "U": 238.02891
}

atomic_symbols = {i: sym for i, sym in enumerate(atomic_masses.keys(), start=1)}

def is_float(s):
    try:
        float(s)
        return True
    except:
        return False


def extract_data_from_log(filename):
    with open(filename, 'r') as f:
        lines = f.readlines()

    atoms = []
    mol_charge = None
    esp_charges = []
    nbo_charges = []
    gaussian_dipole = []

    # --- Extract total charge ---
    for line in lines:
        if 'Charge =' in line and 'Multiplicity =' in line:
            mol_charge = int(float(line.split()[2]))

    # --- Orientation fallback ---
    orientations = ["Standard orientation:", "Z-Matrix orientation:", "Input orientation:"]
    found = None
    block = []

    for orient in orientations:
        blocks = []
        for i, line in enumerate(lines):
            if orient in line:
                tmp = []
                for j in range(i+5, len(lines)):
                    if '-----' in lines[j]:
                        break
                    tmp.append(lines[j].split())
                blocks.append(tmp)
        if blocks:
            block = blocks[-1]
            found = orient
            break

    if found:
        print(f'✅ Using "{found}" for coordinates.')
    else:
        print("❌ No orientation block found. Exiting.")
        sys.exit(1)

    atoms = [(atomic_symbols[int(a[1])], float(a[3]), float(a[4]), float(a[5])) for a in block]

    # --- ESP charges ---
    esp_section = False
    for line in lines:
        if 'ESP charges:' in line:
            esp_section = True
            continue
        if esp_section:
            if not line.strip() or 'Sum of ESP charges' in line:
                break
            parts = line.split()
            if len(parts) >= 3 and is_float(parts[2]):
                esp_charges.append(float(parts[2]))

    # --- NBO charges ---
    nbo_section = False
    for line in lines:
        if 'Summary of Natural Population Analysis:' in line:
            nbo_section = True
        elif nbo_section and '===' in line:
            break
        elif nbo_section and '----' not in line:
            parts = line.split()
            if len(parts) >= 3 and is_float(parts[2]):
                nbo_charges.append(float(parts[2]))

    # --- Gaussian dipole ---
    for i, line in enumerate(lines):
        if 'Dipole moment (field-independent basis, Debye):' in line:
            for j in range(i+1, i+5):
                if 'X=' in lines[j]:
                    parts = lines[j].split()
                    gaussian_dipole = [float(parts[1]), float(parts[3]), float(parts[5])]
                    break

    return atoms, mol_charge, gaussian_dipole, esp_charges, nbo_charges


def compute_center_of_mass(atoms):
    total_mass = sum(atomic_masses[a[0]] for a in atoms)
    weighted_sum = sum(atomic_masses[a[0]] * np.array(a[1:]) for a in atoms)
    return weighted_sum / total_mass


def compute_dipole(atoms, charges, origin):
    dipole = np.zeros(3)
    for i, (_, x, y, z) in enumerate(atoms):
        dipole += charges[i] * (np.array([x, y, z]) - origin)
    return dipole * 4.80320427


def write_xyz_with_dipole(filename, atoms, origin, dipole_vec, label):
    """Appends dummy atoms X (origin) and Bq (dipole head, 6Å away) and saves an XYZ."""
    dipole_dir = dipole_vec / np.linalg.norm(dipole_vec)
    Bq = origin + 6.0 * dipole_dir  # 6Å arrow length. You can change this value as per need
    X = origin

    xyz_lines = [f"{len(atoms) + 2}\n{label} dipole visualization (X,Bq added)\n"]
    for sym, x, y, z in atoms:
        xyz_lines.append(f"{sym:2s} {x:12.6f} {y:12.6f} {z:12.6f}\n")
    xyz_lines.append(f"X {X[0]:12.6f} {X[1]:12.6f} {X[2]:12.6f}\n")
    xyz_lines.append(f"Bq {Bq[0]:12.6f} {Bq[1]:12.6f} {Bq[2]:12.6f}\n")

    outname = os.path.splitext(filename)[0] + f"_{label}_dipole.xyz"
    with open(outname, 'w') as f:
        f.writelines(xyz_lines)
    print(f"💾 Saved dipole visualization: {outname}")


# === MAIN ===
if __name__ == "__main__":
    if len(sys.argv) < 2:
        print("Usage: python3 script.py file.log [-atom N] [-origin X Y Z]")
        sys.exit(1)

    filename = sys.argv[1]
    atom_origin = None
    manual_origin = None

    if "-atom" in sys.argv:
        try:
            atom_origin = int(sys.argv[sys.argv.index("-atom") + 1])
        except:
            print("⚠️ Invalid atom number after -atom.")
            sys.exit(1)

    if "-origin" in sys.argv:
        try:
            i = sys.argv.index("-origin")
            manual_origin = np.array([float(sys.argv[i+1]), float(sys.argv[i+2]), float(sys.argv[i+3])])
        except:
            print("⚠️ Invalid coordinates after -origin (need X Y Z).")
            sys.exit(1)

    atoms, mol_charge, g_dipole, esp, nbo = extract_data_from_log(filename)

    # Determine origin
    if manual_origin is not None:
        origin = manual_origin
        print(f"✅ Using manual origin: {origin}")
    elif atom_origin is not None:
        if 1 <= atom_origin <= len(atoms):
            atom = atoms[atom_origin - 1]
            origin = np.array(atom[1:])
            print(f"✅ Using Atom #{atom_origin} ({atom[0]}) as origin: {origin}")
        else:
            print("⚠️ Atom index out of range.")
            sys.exit(1)
    else:
        origin = compute_center_of_mass(atoms)
        print(f"✅ Using center of mass as origin: {origin}")

    # Compute and save dipoles
    if esp:
        mu_esp = compute_dipole(atoms, esp, origin)
        print(f"ESP Dipole (wrt origin): {mu_esp}  |μ|={np.linalg.norm(mu_esp):.4f}")
        write_xyz_with_dipole(filename, atoms, origin, mu_esp, "ESP")
    else:
        print("⚠️ ESP charges not found.")

    if nbo:
        mu_nbo = compute_dipole(atoms, nbo, origin)
        print(f"NBO Dipole (wrt origin): {mu_nbo}  |μ|={np.linalg.norm(mu_nbo):.4f}")
        write_xyz_with_dipole(filename, atoms, origin, mu_nbo, "NBO")
    else:
        print("⚠️ NBO charges not found.")