#!/Users/lyq/anaconda3/envs/chem/bin/python3

##-------------------------------------##--------------------------------------##
# This is a script which cooparate with "solv_g16.sh"                           #             
# This scipt can extract the energy from xx_solv.log and xx_dry.log, then       #             
# output the Solvation Gibbs Free Energy in kcal/mol (E_solv-E_dry)*627.51+1.87 #                                                                                         
#                                                                               #             
# Y.Q. Li                                  Last update: Jan 18 2025             #             
##-------------------------------------##--------------------------------------##

# HITACHI MAKO IS MY WAIFU! #
import glob,os

def get_last_scf_energy(logfile):
    last_energy = None
    with open(logfile, 'r') as f:
        for line in f:
            if "SCF Done" in line:
                parts = line.split()
                try:
                    energy = float(parts[4])
                    last_energy = energy
                except ValueError:
                    pass
    return last_energy

def main():
    log_files = glob.glob("*.log")
    energies_dict = {}
    for logfile in log_files:
        if logfile.endswith("_dry.log"):
            prefix = logfile[:-8]  
            scf_energy = get_last_scf_energy(logfile)
            energies_dict.setdefault(prefix, {})["dry"] = scf_energy
        elif logfile.endswith("_solv.log"):
            prefix = logfile[:-9]  
            scf_energy = get_last_scf_energy(logfile)
            energies_dict.setdefault(prefix, {})["solv"] = scf_energy
        else:
            continue
    
    for prefix, energy_pair in energies_dict.items():
        e_dry = energy_pair.get("dry", None)
        e_solv = energy_pair.get("solv", None)
        if e_dry is not None and e_solv is not None:
            delta_e = e_solv - e_dry
            delta_g_solv = delta_e * 627.51 + 1.89
            print(f"{prefix}: {delta_g_solv:.4f}")
        else:
            if e_dry is None:
                print(f"{prefix}: PROGRAM ABORTED:No vacuum energy found!")
            if e_solv is None:
                print(f"{prefix}: PROGRAM ABORTED:No solvated energy found!")

if __name__ == "__main__":
    main()