read_resubmit.py

#!/usr/bin/env python3
# -*- coding: utf-8 -*-

"""
Script that reads the output from many Gaussian log files, and prepares a new input file for another software keeping
geometry, charge, multiplicity, etc.
Last Update 2024-03-11 by Emmanuel Nicolas
email emmanuel.nicolas -at- cea.fr
Requires Python3 to be installed.
"""

import argparse
import sys
import cclib
from cclib.parser.utils import PeriodicTable
from ccinput.wrapper import gen_input
from pathlib import Path


def main():
    """
        Main function.
        Retrieves names of files to extract, eventual parameters, starts
        extraction, exports the results in txt files.
    """
    # Get parameters from command line
    runvalues = get_options()
    list_of_raw_files = runvalues['files']
    # Parse all files
    structures = parse_all_files(list_of_raw_files)
    # Use structure data and create the molecules
    write_input_files(structures, runvalues)


def get_options():
    """
        Check command line options and accordingly set computation parameters
    """
    parser = argparse.ArgumentParser(description=help_description(),
                                     epilog=help_epilog())
    parser.formatter_class = argparse.RawDescriptionHelpFormatter
    parser.add_argument('--files', type=str, nargs='+',
                        help='The Gaussian log files from which to extract the relevant energies.')
    parser.add_argument('--software', type=str, default="gaussian", help="Software used for the input")
    parser.add_argument('--type', type=str, help="Type of calculation")
    parser.add_argument('--method', type=str, help="Functional for DFT")
    parser.add_argument('--basis_set', type=str, help="The basis set")
    parser.add_argument('--density_fitting', type=str, help="density_fitting")
    parser.add_argument('--solvent', type=str, help="Solvent")
    parser.add_argument('--solvation_model', type=str, help="solvation model")

    try:
        args = parser.parse_args()
    except argparse.ArgumentError as error:
        print(str(error))  # Print something like "option -a not recognized"
        sys.exit(2)

    runvalues = vars(args)
    # Get values from parser
    runvalues['files'] = [Path(x).resolve() for x in args.files]
    return runvalues


def extract_data(file):
    energies = dict()

    data = cclib.io.ccread(file)
    energies["xyz"] = xyz_geometry(data.atomcoords[-1], data.atomnos.tolist())
    energies["charge"] = data.charge
    energies["multiplicity"] = data.mult

    return energies


def xyz_geometry(coordinates, elements):
    """Returns geometry in XYZ format"""
    periodic_table = PeriodicTable()
    xyz = [
        " ".join(
            [periodic_table.element[elements[i]].ljust(5)]
            + ["{:.9f}".format(s).rjust(25) for s in atom]
        )
        for i, atom in enumerate(coordinates)
    ]
    return xyz


def parse_all_files(files):
    """
        Parse all files in list of files
    """
    structures = list()

    for file in files:
        data = extract_data(file)
        data["name"] = file
        structures.append(data)

    return structures


def write_input_files(structures, runvalues):
    """
        For each structure, write the input file using the ccinput engine as a library
    """
    for molecule in structures:
        inp = gen_input(software=runvalues["software"],
                        type=runvalues["type"],
                        xyz="\n".join(molecule["xyz"]),
                        name=molecule["name"].stem,
                        aux_name=molecule["name"].stem,
                        method=runvalues["method"],
                        charge=molecule["charge"],
                        multiplicity=molecule["multiplicity"],
                        basis_set=runvalues["basis_set"],
                        density_fitting=runvalues["density_fitting"],
                        solvent=runvalues["solvent"],
                        solvation_model=runvalues["solvation_model"])
        input_file = molecule["name"].with_suffix(".com")
        with open(input_file, mode="w") as file:
            file.write(inp)


def help_description():
    """
        Returns description of program for help message
    """
    return ""


def help_epilog():
    """
        Returns additional help message
    """
    return ""


if __name__ == '__main__':
    main()