How to estimate binding energy with Dreiding Forcefield using frozen pose method

[ScottLiao920]

Hi, thanks for your amazing work!

Just curious about the actual code for binding energy estimation with Dreiding Forcefield using frozen pose method.
Suppose I already got one docked complex of zeolite and OSDA stored in .cif file, is the following way to estimate binding energy correct?

from gulpy.jobs import OptiJob
from gulpy.structure import GulpCrystal
from gulpy.inputs import Library
from gulpy.parser import PropertyParser, JobParser
from gulpy.structure import GulpCrystal, GulpMolecule
from gulpy.structure.labels import CatlowLabels, DreidingMoleculeLabels, DreidingLabels
from pymatgen.core import Structure, Molecule, Lattice

structure = Structure.from_file(cif_path) # cif for the docked complex
structure_gulpy = GulpCrystal(structure, DreidingLabels())
zeolite_gulpy = GulpCrystal(
    Structure.from_file(zeolite_cif_path),
    DreidingLabels(),
)
osda_gulpy = GulpMolecule.from_smiles(
    Molecule.from_file(osda_conformer_xyz_path).cart_coords,
    osda_smiles,
    DreidingMoleculeLabels(),
)
library = Library.from_file(f"{os.environ["GULP_LIB"]}/dreiding.lib")

# Create optimization job with defaults
job = OptiJob(
    structure_gulpy,
    library,
    keywords=["opti", "conp"],  # Optimize it under constant pressure
    options={"maxcyc": 300},  # Increase max cycles
)

# Write input file
job.write_input("optimization.gin")

# Run optimization
job.run("optimization.gin")

# Parse results
results = job.parse_results("optimization.out")

# repeat for zeolite, OSDA with the same settings

and the estimated binding energy is results['energy']-zeo_results['energy'] - osda_result['energy'] in eV?
I am not sure about whether my optimization job parameters or the library is right or not, please correct me if I'm wrong.

Thanks for your help in advance.