Adding support for Nudged Elastic Band calculations

[TroyNSmith]

@avcopan and I are gearing up to implement NEB calculations through geomeTRIC. We wanted to have your input on what that should look like. The main difference from a standard calculation being that there are multiple input geometries: starting_structure, ending_structure, and possibly guess_saddle_structure.

[coltonbh]

Great questions. I'll check in with @jandestrada in my lab on Tuesday to chat about this. He's been developing some NEB software as well and has built various data structures for it on top of qcdata data structures. He and I will be doing more NEB/reaction path work over the coming months, so we need to come up with something general purpose that can work long term. Will circle back with ideas later this week :)

[avcopan]

Thanks @coltonbh ! Looking forward to it!

[coltonbh]

Just circling back on this, @avcopan and @TroyNSmith... Am teaching this week and at a conference next week. Hoping to get to this over the weekend. Are you currently blocked on this? Gauging time urgency here!

[TroyNSmith]

@coltonbh Sorry for my delayed response! I'll be tied up with graduations through this afternoon and all day tomorrow, so if you're able to give us a sense of direction this weekend, that would be fine! Otherwise we've got some other tasks I can shift my attention back to if you need more time, so please don't feel rushed!

[coltonbh]

Ok great! @TroyNSmith I'll keep you posted on what I can do so as to not leave you hanging! A lot happening for graduation over here and finishing up in my lab... Excited to be full time on this software stack soon :)

[coltonbh]

Ok I have a proposal, I'd be curious for your thoughts @TroyNSmith and @avcopan:

Option 1

class ProgramInput(...):
    structure: Structure | None = None
    structures: dict[str, Structure] | None = None

Single-point, gradient, hessian, optimization, conformer search, etc. continue to use:

ProgramInput(structure=water, ...)

NEB-like jobs use:

ProgramInput(
    calctype="neb",
    structures={
        "reactant": reactant,
        "product": product,
      },
      ...
  )

This looks nicest for NEB style calculations, but we loose type safety for .structure always being present instead of Structure | None.

Option 2

class ProgramInput(...):
    structure: Structure
    structures: dict[str, Structure] = {}

Then define the invariant:

• .structure is always the primary/default/start structure.
• .structures contains additional named structures needed by specialized calculations.
• For single-structure jobs, .structures == {}.
• For NEB:
  • .structure is the reactant/start structure.
  • .structures["product"] is the product/end structure.

This preserves type safety and backwards compatibility. It also keeps the common case clean.

The tradeoff is that NEB becomes slightly asymmetric: reactant lives at .structure, product lives in .structures["product"]. I think that is acceptable because many path/trajectory calculations naturally have a “starting” structure, and .structure already means “the structure this calculation starts from” throughout the codebase. Then we document: input.structure means “the primary/start/reference structure.” input.structures means “additional named structures required by this calculation.”

For something like constrained docking or embedding it would look like this:

ProgramInput(
    calctype="optimization",
    structure=guest,
    structures={"host": host},
    model=...,
)

For surface chemistry:

ProgramInput(  
    calctype="optimization",
    structure=adsorbate_on_surface,
    structures={
        "clean_surface": surface,
        "isolated_adsorbate": adsorbate,
    },
    model=...,
)

I prefer option 2 so that we don't loose type safety for all single-structure calculations. Thoughts? Other ideas?

Once we decide on something I'll get it implemented and released in a day so you guys can proceed.

[avcopan]

@coltonbh Thanks for getting back to us! I prefer option 2 as well, and I think this should be fairly intuitive to users as it is also the way it works in QC programs, at least in Orca.

[avcopan]

So from our end you are good to go ahead with the implementation. This will be great to have!

[coltonbh]

Great! Will implemented and release tomorrow :)

[TroyNSmith]

@coltonbh Thank you so much for the time and consideration you've put into helping this move forward! I hope that your graduation and all related activities went (/ are going) smoothly!

I know that I'm late here but I agree with @avcopan on option 2 being the best option. I'm wondering if this layout can (or should) be mirrored in the program output so that we retrieve, from a free-end NEB for example:

...

prog_inp = ProgramInput(
    calctype="neb",
    structure=reactant
    structures={"product": product, "ts_guess": ts_guess},
    model=...,
    keywords={"free_end": True},
)

prog_out = compute("geometric", prog_inp)

final_reactant = prog_out.data.structures["reactant"]
final_product = prog_out.data.structures["product"]
final_transition = prog_out.data.structures["ts_guess"]

...

Perhaps this question comes at a higher-level of implementation or the ProgramOutput is already flexible in this way--I'm mostly placing this here so that questions guiding support for NEB are contained in one public location. :)

[avcopan]

I think a regular sequence of geometries should be fine for the output. The NEB doesn't itself identify the TS, and the end points are easy to identify as the first and last items in the sequence.

[coltonbh]

@avcopan @TroyNSmith

Multi-structure inputs now implemented in qcdata #103 and released as qcdata v0.17.2.

For the Data object that will hold the outputs, why don't you guys look at neb-dynamics, a NEB package my lab has worked on that has build data structures based on qcdata for some ideas? You guys are doing the NEBs so I think you'll know more precisely than me how to build a good data structure for holding the outputs. Can you guys design that and post a proposal here? As an added benefit it will probably spill over in the neb-dynamics package, which will be nice. Once we agree to something, if you guys can implement it and make sure it's working for your cases well, that would be great!

I'm happy to offer perspective and be quick on design reviews :) If you need more support than this, please let me know!