#!/usr/bin/env python
# coding: utf-8

# ## Initialization

from scm.plams import Settings, Molecule, AMSJob, Units, view, plot_image_grid
from scm.plams.recipes.adffragment import ADFFragmentJob


# ## Define the molecules
# For convenience we define here two molecules, normally you would read them from xyz files


def get_molecule(input_string):
    job = AMSJob.from_input(input_string)
    return job.molecule[""]


mol1 = get_molecule(
    """
System
    Atoms
        C      -0.75086900       1.37782400      -2.43303700
        C      -0.05392100       2.51281000      -2.41769100
        H      -1.78964800       1.33942600      -2.09651100
        H      -0.30849400       0.43896500      -2.76734700
        H      -0.49177100       3.45043100      -2.06789100
        H       0.98633900       2.54913500      -2.74329400
    End
End
"""
)


mol2 = get_molecule(
    """
System
    Atoms
        C       0.14667300      -0.21503500       0.40053800
        C       1.45297400      -0.07836900       0.12424400
        C       2.23119700       1.15868100       0.12912100
        C       1.78331500       2.39701500       0.38779700
        H      -0.48348000       0.63110600       0.67664100
        H      -0.33261900      -1.19332100       0.35411600
        H       2.01546300      -0.97840100      -0.14506700
        H       3.29046200       1.03872500      -0.12139700
        H       2.45728900       3.25301000       0.35150400
        H       0.74193400       2.60120700       0.64028800
    End
End
"""
)

# Molecules could also be read from xyz files
# mol1 = Molecule("ethene.xyz")
# mol2 = Molecule("butadiene.xyz")

plot_image_grid(
    {"1": view(mol1, guess_bonds=True), "2": view(mol2, guess_bonds=True)}, show_labels=False, save_path="picture1.png"
)
# ## Setup and run the job

common = Settings()  # common settings for all 3 jobs
common.input.ams.Task = "SinglePoint"
common.input.adf.basis.type = "DZP"
common.input.adf.xc.gga = "PBE"
common.input.adf.symmetry = "NOSYM"

full = Settings()  # additional settings for full system calculation
full.input.adf.etsnocv  # empty block
full.input.adf.print = "etslowdin"

job = ADFFragmentJob(fragment1=mol1, fragment2=mol2, settings=common, full_settings=full)
results = job.run()


# ## Energy Decomposition Results

# Energy terms contributing to the bond energy (with respect to the fragments):

import pandas as pd

bond_energy = results.get_energy(unit="eV")
decom = results.get_energy_decomposition(unit="eV")

df = pd.DataFrame(
    [{"Term": t, "Energy [eV]": e} for t, e in list(decom.items()) + [("Total Bond Energy", bond_energy)]]
)
print(df)