[FEATURE] - Add optimizations with Sella/ASE

[avcopan]

I like the API for qcop and would love to use it as a common interface for energy, hessian, and structure calculations with Sella/ASE, similar to the current geomeTRIC/QCEngine calculators in qcop. The motivation for adding the former is two-fold:

1. Sella seems to have more robust support for saddle-point/TS optimizations
2. ASE provides a QCEngine-style wrapper around another set of electronic structure codes, in addition to those already present in qcop

In particular, my workflow involves the following list of routines, which might be performed through Sella/ASE:

• ASE: Single-point energy calculations using an ASE calculator (ideally, preserving the option to pass in a custom calculator)
• ASE: Gradient (force) calculations using an ASE calculator
• ASE: Hessian (harmonic vibration) calculations using an ASE calculator
• Sella: Optimization of a minimum-energy (equilibrium) structure using internal coordinates
• Sella: Optimization of a saddle-point (TS) structure using internal coordinates
• Sella: Optimization of a minimum energy structure, subject to one or more constraints of bond distance, angle, etc., also in internal coordinates (for generating relaxed scans)
• Sella: IRC calculation in internal coordinates (lower priority than the other items, but would be nice to have)

I don't need any further tools from ASE, beyond these basic ones. These building blocks would allow me to write my own routines for conformer sampling, relaxed scanning, and transition state finding using qcop.

[coltonbh]

Thanks for the detailed write up! Super helpful.

A few clarifications:

1. Is it important that your single point/gradient/hessian/frequency calculations run via ASE? Or just that they get run with the requested package (e.g., TeraChem, QChem, etc.) using a qcop interface? Are the packages you want to use for single point calcs already covered in qcop or are specific packages/calctypes missing that you need? If some are missing, which ones? Based on what you wrote here, my impression is they are all covered, but I want to check.

Sella does seem pretty slick! My initial impression for how I'd implement this is provide the usual qcop interface for optimizations and transition state calcs for Sella and then add a qcop ASE calculator that runs the single point calculations for powering the optimization routine using qcop so that we collect all of that data back in a nice, structured way. This is why I ask about the ASE calculations you list--is it important that we wrap up specific features of ASE or is it just important that these calculations get done (via a usual qcop calculation routine) both as individual calculations and as part of Sella's optimization routines?

Thanks for the detail and time to write this up. As they say--measure twice, cut once :) So I try to get clarity on exactly what needs to be built before getting started :sparkles:

[avcopan]

No, it isn't important that those things run via ASE. Any way I can get them is fine, as long as the results are consistent (i.e. running a single point on the optimized geometry gives me the same answer).

[avcopan]

Regarding electronic structure back-ends, the main packages I plan to use in the near term I believe are all available in qcop. Off the top of my head: PySCF, Psi4, Molpro, Orca, XTB. I would be fine with simply adding packages to qcop as I encounter the need. I've become fairly decent with writing my own parsers over the years.

If it were possible, though, it would be cool if the interface could be done in a generic way, so that all of the ASE back-ends are automatically available through qcop. But I gather that that is not possible because the ASE calculators are missing information that qcop collects. Correct?

Another "nice to have" feature from ASE is the ability to externally define a "calculator" and pass it into the program. I really like this feature, because of the flexibility it would provide to users. For example, in my PhD I implemented my own electronic structure methods. It would have been great to be able to simply define some sort of "calculator" object and instantly be able to plug into any routines that were built on top of qcop. This could be done independently of ASE, and even if it took a little bit more work to provide all of the appropriate information that qcop needs, I think it could be worth it.

But again, those are separate features. Simply being able to run Sella through qcop would be a huge selling point for me.

[coltonbh]

Yes. I could add an ASEAdapter, like the QCEngineAdapter I have in qcop, that generically allows us to use everything in ASE. That's probably a good way to go.

You can already add new Adapters dynamically, like you described with ASE, by following the example code below :) If this is missing something you had in mind, please let me know. This allows you to create new methods in your own code and "register" them with qcop without ever needing to change the source code :man_technologist:

[coltonbh]

Here's a quick version; I should add some better docs for this. This also makes me want to move the update_func stuff up a level so end users don't have to consider it.

All you have to do is create a new Adapter that inherits from ProgramAdapter. This may look a little fancy with the typing Generics (those just give you full type safety), but all you need to do is:

1. Fill in the supported_calctypes with the calculation types your Adapter will support.
2. Give the Adapter the program name.
3. Implement program_version.
4. Implement compute_results.

With that, qcop knows about your new program and is ready to use it :)

from typing import Callable, Optional, Tuple

from qcio import CalcType, ProgramInput, SinglePointResults, Structure

from qcop import compute, exceptions
from qcop.adapters.base import ProgramAdapter

class MyAdapter(ProgramAdapter[ProgramInput, SinglePointResults]):
    # Fill in supported_calctypes and program
    supported_calctypes = [CalcType.energy]
    program = "myprog"

    def program_version(self, stdout: Optional[str]) -> str:
        """Modify to return the actual version of the program.

        It can find this value from the command line, a Python package, or by parsing
        the stdout from the program.
        """
        return "v0.0.1"

    def compute_results(
        self,
        inp_obj: ProgramInput,
        update_func: Optional[Callable] = None,
        update_interval: Optional[float] = None,
        **kwargs
    ) -> Tuple[SinglePointResults, str]:
        """Perform the calculation with external Program. Return the results and stdout."""
        # Perform the calculation and return the results and stdout
        return SinglePointResults(energy=123.456), "stdout stuff"

# Now your new Adapter is ready to use inside of qcop!
prog_input = ProgramInput(
    calctype="energy",
    structure=Structure(symbols=["H", "O"], geometry=[[0, 0, 0], [0, 0, 1]]),
    model={"method": "HF", "basis": "sto-3g"},
)

try:
    po = compute("myprog", prog_input)
except exceptions.QCOPBaseError as e:
    po = e.program_output

[avcopan]

Wow, that sounds great! If it all works, this would be my dream! I have been looking for a package that does these things literally for years.

As soon as you have an interface for Sella optimizations, I can start prototyping code on top of qcop and submitting issues or (for simple things) submitting PRs.

[coltonbh]

@avcopan sorry I've been slow on this! Added a new feature this week that I think you'll find really valuable--the ability to visualize all your results in a Jupyter notebook with a single command! Details here. I needed this for viewing some of my own results this week so I prioritized it. Sella/ASE up next...

[avcopan]

No worries! I am in for a busy month, so it isn't a huge rush from my end. The visualization looks nice! Are you using py3Dmol or something else?

[coltonbh]

Yes, py3dmol for the molecular viewers. Would you suggest trying out anything else? It works pretty well except it doesn't handle a large number of viewers at the same time super well. The developer is adding some handy fixes though which make that a little smoother. And then of course matplotlib for the plots and a little custom HTML/CSS to create the tables.

Ok good to hear it's not holding you up! Will get to it in the next few weeks then :)

[avcopan]

No, I like py3dmol! It's what I use as well.