Mod:Hunt Research Group/center of mass

Code to Calculate the Center of Mass (CoM) and produce xyz file with CoM

What the program does?

1. The script can read an XYZ or Gausslain log file and extract its coordinates

2. Then it calculates the center of mass (CoM) from those coordinates

3. It then saves the file_with_COM.xyz to the current directory with a dummy atom (X) with CoM corordinates

Running the python Script

• Running the script is very simple. Simply save the script as center-of-mass_calculation.py and run it as follows (considering your molecule file is water.log)

python3 center-of-mass_calculation.py water.log

• It will produce a file named water_with_CoM.xyz that will have a dummy atom X for showing the center of mass.

Below is the Script

#!/usr/bin/env python3
"""
Script to calculate the center of mass (CoM) of a molecule from an XYZ/log file,
and write a new XYZ file with the CoM added as a pseudo-atom 'X'.
Script by Dr. Muhammad Ali Hashmi (16/05/2025)
"""

import sys
import os

# Atomic weights dictionary
symbol_to_mass = {
"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, "Np": 237.0, "Pu": 244.0,
"Am": 243.0, "Cm": 247.0, "Bk": 247.0, "Cf": 251.0, "Es": 252.0, "Fm": 257.0,
"Md": 258.0, "No": 259.0, "Lr": 262.0, "Rf": 267.0, "Db": 268.0, "Sg": 271.0,
"Bh": 272.0, "Hs": 270.0, "Mt": 276.0, "Ds": 281.0, "Rg": 280.0, "Cn": 285.0,
"Nh": 286.0, "Fl": 289.0, "Mc": 290.0, "Lv": 293.0, "Ts": 294.0, "Og": 294.0}

# Atomic number to symbol mapping
number_to_symbol = {i: symbol for i, symbol in enumerate(symbol_to_mass.keys(), start=1)}

# A function to read the coordinates from an xyz file
def read_xyz_coordinates(filepath):
    coords = []
    with open(filepath, 'r') as f:
        lines = f.readlines()[2:]
        for line in lines:
            parts = line.split()
            if len(parts) >= 4:
                symbol = parts[0]
                x, y, z = map(float, parts[1:4])
                coords.append([symbol, x, y, z])
    return coords

# A function to read the coordinates from a gaussian log file
def read_log_coordinates(filepath):
    coords = []
    with open(filepath, 'r', encoding='utf-8', errors='ignore') as f:
        lines = f.readlines()
    for i, line in enumerate(lines):
        if "Standard orientation:" in line:
            coords.clear()
            start = i + 5
            while start < len(lines):
                if "-----" in lines[start]:
                    break
                parts = lines[start].split()
                atomic_num = int(parts[1])
                symbol = number_to_symbol.get(atomic_num)
                x, y, z = map(float, parts[3:6])
                coords.append([symbol, x, y, z])
                start += 1
    return coords

# A function to calculate the center of mass
def calculate_center_of_mass(coords):
    x_cm = y_cm = z_cm = total_mass = 0.0
    for symbol, x, y, z in coords:
        mass = symbol_to_mass.get(symbol)
        if mass is None:
            raise ValueError(f"Unknown atomic symbol: {symbol}")
        x_cm += mass * x
        y_cm += mass * y
        z_cm += mass * z
        total_mass += mass
    return x_cm / total_mass, y_cm / total_mass, z_cm / total_mass

# A function to write the molecule with CoM coordinates as a dummy atom X
def write_xyz_with_CoM(filename, coords, CoM):
    output_file = os.path.splitext(filename)[0] + "_with_CoM.xyz"
    with open(output_file, 'w') as f:
        f.write(f"{len(coords)+1}\n")
        f.write(f"Center of Mass added as dummy atom X\n")
        for atom in coords:
            f.write(f"{atom[0]:<2} {atom[1]:>12.6f} {atom[2]:>12.6f} {atom[3]:>12.6f}\n")
        f.write(f"X  {CoM[0]:>12.6f} {CoM[1]:>12.6f} {CoM[2]:>12.6f}\n")
    print(f"Written: {output_file}")

# === Main ===
if __name__ == "__main__":
    if len(sys.argv) < 2:
        print("Usage: python center_of_mass.py <file.xyz or .log>")
        sys.exit(1)

    filename = sys.argv[1]

    if not os.path.isfile(filename):
        print(f"Error: File '{filename}' not found.")
        sys.exit(1)

    try:
        if filename.endswith(".xyz"):
            coords = read_xyz_coordinates(filename)
        elif filename.endswith(".log"):
            coords = read_log_coordinates(filename)
        else:
            raise ValueError("Unsupported file format. Use .xyz or .log only.")

        CoM = calculate_center_of_mass(coords)
        print(f"{filename}: Center of Mass = ({CoM[0]:.3f}, {CoM[1]:.3f}, {CoM[2]:.3f})")

        write_xyz_with_CoM(filename, coords, CoM)

    except Exception as e:
        print(f"Error: {e}")
        sys.exit(1)