Just pointing out that this PR doesn't remove C-side optking code yet. Is that coming?

Just pointing out that this PR doesn't remove C-side optking code yet. Is that coming?

Just an oversight. deleted but not comitted.

Just pointing out that this PR doesn't remove C-side optking code yet. Is that coming?

Just an oversight. deleted but not comitted.

The Optking is dead. Long live the Pyoptking.

The test that is currently failing is a test in gcp/pbeh3c/. This is due to an optking side issue where the CustomHelper class being used by optking is not accepting a psi4.core.Molecule the type checking was looking for qcdb.Molecule. As a backup optking defaulted to psi4's active molecule.

This is the call.
E = optimize('pbeh3c/def2-msvp', molecule=unopethene)

optking takes the molecule here in the driver and uses a default fallback instead:
opt_object = optking.opt_helper.CustomHelper(molecule, params=optimizer_params)

The type check will get changed in optking to include core.Molecule. An alternative question this raises for me is whether the active_molecule should get updated at some point in the optimization. Updating the active molecule patches the issue but is that desired?

The test that is currently failing is a test in gcp/pbeh3c/. This is due to an optking side issue where the CustomHelper class being used by optking is not accepting a psi4.core.Molecule the type checking was looking for qcdb.Molecule. As a backup optking defaulted to psi4's active molecule.

This is the call. E = optimize('pbeh3c/def2-msvp', molecule=unopethene)

optking takes the molecule here in the driver and uses a default fallback instead: opt_object = optking.opt_helper.CustomHelper(molecule, params=optimizer_params)

The type check will get changed in optking to include core.Molecule. An alternative question this raises for me is whether the active_molecule should get updated at some point in the optimization. Updating the active molecule patches the issue but is that desired?

My vote would be "no"; it's too much like changing the input. However, I do think that some users would expect the final active molecule to be updated for them. And workflows inside the python input may seem more intuitive that way.

An alternative question this raises for me is whether the active_molecule should get updated at some point in the optimization. Updating the active molecule patches the issue but is that desired?

My vote would be "no"; it's too much like changing the input. However, I do think that some users would expect the final active molecule to be updated for them. And workflows inside the python input may seem more intuitive that way.

I may not be following this right. I'd say the state of the psi4 active mol during an optimization is arbitrary -- whatever works for you. In cpp-optking, I think the communication was through legacymolecule anyways. The molecule optking is acting on should be updated by the time control returns to the user at the end of the opt. I thought this was already happening through https://github.com/psi4/psi4/pull/2727/files#diff-acf663ccea13592c4c656cf89c7b62e6f5bd3b2e8b4f12ba354129bd39d096f8R1296-R1297 . That's consistent with cpp-optking, and I think that must be happening b/c many of the tests check NRE before and after opt.

An alternative question this raises for me is whether the active_molecule should get updated at some point in the optimization. Updating the active molecule patches the issue but is that desired?

My vote would be "no"; it's too much like changing the input. However, I do think that some users would expect the final active molecule to be updated for them. And workflows inside the python input may seem more intuitive that way.

I may not be following this right. I'd say the state of the psi4 active mol during an optimization is arbitrary -- whatever works for you. In cpp-optking, I think the communication was through legacymolecule anyways. The molecule optking is acting on should be updated by the time control returns to the user at the end of the opt. I thought this was already happening through https://github.com/psi4/psi4/pull/2727/files#diff-acf663ccea13592c4c656cf89c7b62e6f5bd3b2e8b4f12ba354129bd39d096f8R1296-R1297 . That's consistent with cpp-optking, and I think that must be happening b/c many of the tests check NRE before and after opt.

OK, you changed my hasty mind on that. What concerns me is that the user may, in some instances, not realize that the default active molecule is the one that optking is acting on and changing. But the upside convenience wins, I agree.

OK, you changed my hasty mind on that. What concerns me is that the user may, in some instances, not realize that the default active molecule is the one that optking is acting on and changing. But the upside convenience wins, I agree.

If changing the mol over the course of an opt was new behavior, I think it'd be very much up for debate. I know FAE is not keen on background/globals changes. And they are weird in a Jupyter notebook. But for now, I'm for minimal surprise upon cpp-optking -> py-optking.

Please ping me for review once Lori approves. She knows this part of the code better than I do.

In cases where reference values are not matched perfectly (but the test should pass) is it better to loosen the comparison or update the reference value. I assume updating the reference value is the way to go but I want to have some record of asking before I start slightly changing reference values.

Secondarily, should the reference values be updated in general at some point so that users don't stumble across an instance where the value is slightly off and wonder why?

Explanation:
I've expanded the number of tests I'm running since I started cleaning up core.cc. In the opt specific tests the convergence is usually tight enough that the nuclear repulsion energies match the reference values just fine for both optimizers. There are some tests like cc1-3 that are failing due to being just above threshold. atol = 0.001 the difference is ~ 0.0017. These tests are using the default qchem convergence criteria ~ 3e-4 max_force. Both optimizers finish well below the criteria and geometries match to 1e-4 Angstroms and 0.001 degrees. Geometries are the same. There are around 5 or 6 tests failing like this.

Would y'all like this PR to include updated output.ref files for the tests or would a separate test updating PR be better? I don't see anything in the" adding tests documentation" about reference log files. I can include those as well if desired. It will just greatly increase the number of lines changed in this PR.

Would y'all like this PR to include updated output.ref files for the tests or would a separate test updating PR be better? I don't see anything in the" adding tests documentation" about reference log files. I can include those as well if desired. It will just greatly increase the number of lines changed in this PR.

Please update .ref files. Probably not worth adding .log files, but that's a @loriab question.

There are some tests like cc1-3 that are failing due to being just above threshold. atol = 0.001 the difference is ~ 0.0017. These tests are using the default qchem convergence criteria ~ 3e-4 max_force. Both optimizers finish well below the criteria and geometries match to 1e-4 Angstroms and 0.001 degrees.

Could you elaborate on why tests are failing at all, and what numbers are differing? Are these Cartesian coordinates?

In cases where reference values are not matched perfectly (but the test should pass) is it better to loosen the comparison or update the reference value. I assume updating the reference value is the way to go but I want to have some record of asking before I start slightly changing reference values.

Secondarily, should the reference values be updated in general at some point so that users don't stumble across an instance where the value is slightly off and wonder why?

General guidance to for ref values to be from a tightly converged/optimized calc, then loosen the comparison check to accommodate the default/existing conv crit (https://psicode.org/psi4manual/master/add_tests.html#test-contents). Tests checking opt status at a certain cycle exempt of course. That's the principle, but do feel free to change as you see fit --- the reference values (agreed, preferred thing to change if the ref is the culprit) or the comparison crit (if it's the optimizer behavior that's the instigator).

Explanation:
I've expanded the number of tests I'm running since I started cleaning up core.cc. In the opt specific tests the convergence is usually tight enough that the nuclear repulsion energies match the reference values just fine for both optimizers. There are some tests like cc1-3 that are failing due to being just above threshold. atol = 0.001 the difference is ~ 0.0017. These tests are using the default qchem convergence criteria ~ 3e-4 max_force. Both optimizers finish well below the criteria and geometries match to 1e-4 Angstroms and 0.001 degrees. Geometries are the same. There are around 5 or 6 tests failing like this.

Thanks for the explanation. I'd view reference NRE values as less venerable. For one thing, only those that caused trouble were even updated when physical constants changed, iirc.

Would y'all like this PR to include updated output.ref files for the tests or would a separate test updating PR be better? I don't see anything in the" adding tests documentation" about reference log files. I can include those as well if desired. It will just greatly increase the number of lines changed in this PR.

Separate, please. Like updating samples/, better to keep the not-for-visual-inspection changes aside. Update: I see Jonathon thinks differently. Including is ok with me, now that GH allows files to be folded up rather than scrolled through. Is .log where all the optking detailed output goes now? I guess we ought to start collecting them. output.log, perhaps. But this can also be deferred to a grand regenerate-the-refs script and PR.

I'm perfectly happy to leave .ref to a separate PR, so long as @AlexHeide agrees to take responsibility for getting that in by PsiCon at the latest. 🙂

Could you elaborate on why tests are failing at all, and what numbers are differing? Are these Cartesian coordinates?

I poorly worded my explanation. I was attempting to say that the geometries are virtually the same but not numerically, they are both well converged. I assumed the discrepancy was just a numeric difference and there have been various tweaks as well over the last few years to the algorithm.

There is 1 algorithmic difference I'd be concerned about, from looking at the two outputs. New optking doesn't consider the most recent step in the hessian updating procedure - cpp-optking did. I've found three pieces of logic that explicitly or implicitly prevent updating with the most recent step in all or specific cases. I'd have to ask @psi-rking if this is a bug or was changed due to some stability concern.

Is .log where all the optking detailed output goes now? I guess we ought to start collecting them. output.log, perhaps. But this can also be deferred to a grand regenerate-the-refs script and PR.

Yes optking's detailed output goes to .log but it isn't strictly speaking necessary for the user to see the detailed logs I would say. However, if the test is failing it might be nice to have a more detailed optimization record to compare against if the optimizer is at fault.

Just updating the reference is fine with me, in that case.

There is 1 algorithmic difference I'd be concerned about, from looking at the two outputs. New optking doesn't consider the most recent step in the hessian updating procedure - cpp-optking did. I've found three pieces of logic that explicitly or implicitly prevent updating with the most recent step in all or specific cases. I'd have to ask @psi-rking if this is a bug or was changed due to some stability concern.

Interesting. This does not ring a bell for me. I take the question to be "given a hessian and the current forces (beyond the first step) should you update the hessian with those forces before using them to calculate the step?" Can you point to the code? IDK, probably best answer is whatever works better in practice and doesn't cause problems.

Actually, it's possible I did this because I learned to avoid hessian updating when displacements are very small (or the geometries are very close). Perhaps I didn't want to update hessian until I knew the step size.

I may not be following this right. I'd say the state of the psi4 active mol during an optimization is arbitrary -- whatever works for you. In cpp-optking, I think the communication was through legacymolecule anyways. The molecule optking is acting on should be updated by the time control returns to the user at the end of the opt. I thought this was already happening through https://github.com/psi4/psi4/pull/2727/files#diff-acf663ccea13592c4c656cf89c7b62e6f5bd3b2e8b4f12ba354129bd39d096f8R1296-R1297 . That's consistent with cpp-optking, and I think that must be happening b/c many of the tests check NRE before and after opt.

I want to make sure that the active molecule behavior is as expected. Whatever molecule the driver uses will be updated. This will be either the active molecule OR the passed molecule. If the molecule is passed the active molecule is not updated in any way. All the asserts pass here.

import math

molecule h2o {
    pubchem:water
}

molecule h2o2 {
    pubchem:hydrogen peroxide
}

# quick comparison. h2o2 is active molecule
h2o2_nre = h2o2.nuclear_repulsion_energy()
active_nre = core.get_active_molecule().nuclear_repulsion_energy()

assert math.isclose(h2o2_nre, active_nre)

# optimize h2o2 (active molecule) expect repulsion energy to match
E = optimize("scf/sto-3g")
h2o2_opt_nre = h2o2.nuclear_repulsion_energy()
active_nre = core.get_active_molecule().nuclear_repulsion_energy()

assert math.isclose(h2o2_opt_nre, active_nre)

# optimize h2o. nuclear repulsion does not match. active molecule is still h2o2
# Currently (next commit will fix) h2o2 would be optimized by this call. (optking side issue)
E = optimize("scf/sto-3g", molecule=h2o)
h2o_opt_nre = h2o.nuclear_repulsion_energy()
active_nre = core.get_active_molecule().nuclear_repulsion_energy()

assert math.isclose(h2o2_opt_nre, active_nre)

The whole wrong molecule being optimized thing is an optking side issue that is fixed on optking/master.

Let me be more explicit about why I find the "multi-fragment optimizations" section confusing.

• The first sentence talks about "the metric for connecting atoms" without explaining what "connecting atoms" means. Even worse, the first sentence is not obviously about dimers.
• It isn't clear to me what a "reference point" signifies. While I can tell whether something is an acceptable reference point, what are these used for? It looks like these are atoms used to define intermolecular internal coordinates.
• I don't know what it means to talk about a linear combination of atoms, or how to interpret [[3], [1], [2]], [[1, 2, 3, 4, 5, 6], [2], [6]] in your first example.

The third point is crucial - I can't follow your examples. I won't insist on a figure, but I do insist on examples I can understand.

looks like typos in options snagging up the docs:

/home/runner/work/psi4/psi4/code/objdir/doc/sphinxman/source/optking.rst:359: ERROR: Undefined substitution referenced: "optking_frag_ref_atoms".
/home/runner/work/psi4/psi4/code/objdir/doc/sphinxman/source/optking.rst:359: ERROR: Undefined substitution referenced: "otking__frag_ref_atoms".