[FEATURE] - Add optimizations with Sella/ASE
Opened this issue ยท 10 comments
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:
- Sella seems to have more robust support for saddle-point/TS optimizations
- 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.
Thanks for the detailed write up! Super helpful.
A few clarifications:
- 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 inqcop
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 โจ
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).
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.
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 ๐จโ๐ป
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:
- Fill in the
supported_calctypes
with the calculation types yourAdapter
will support. - Give the
Adapter
theprogram
name. - Implement
program_version
. - 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
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.
@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...
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?
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 :)
No, I like py3dmol! It's what I use as well.