Using Psi4 as an external optimizer

I would like to use Psi4’s optimizer for geometry optimizations. I have a python (and perl) script which can compute energy and gradients. I couldn’t find much regarding external optimization in the manual. Is it possible to use Psi4’s optimizer while calculating gradients using a different program?

See Geometry Optimisation with arbitrary wavefunction.

Are there special features you need?

Could be easier to use: https://github.com/psi-rking/optking

Thanks a lot for your suggestions. At the moment, I don’t think I need any special features. I will try both methods but it seems like using optking (https://github.com/psi-rking/optking) is a better option, since I just need an optimizer.

I am new to Psi4 and I have been trying your method, but it seems like the function I have written is not returning the gradient in the format that psi4 wants. I am calling the gradient function as follows, psi4.optimize(grad_fun.get_grad(jobname, geom, elem, mcharge, multiplicity), dertype=1, molecule=test)

The function get_grad takes geometry, molecular charge, and multiplicity as inputs and returns a numpy array of the gradient, and a random quantity (charge in my case) as the second argument. I am getting this error, “AttributeError: ‘tuple’ object has no attribute ‘lower’” for line 986 in driver.py

The first argument to optimize needs to be a gradient function, not the result of evaluating a gradient function.

I apologize but I am not sure if I understood what exactly you meant. Do I need to call my gradient function as follows then,
psi4.optimize(opt_func=grad_fun.get_grad(jobname, geom, elem, mcharge, multiplicity), dertype=1, molecule=test)
The above syntax seems to be working. But optimize requires name argument as well. I am not sure what value to give for that since I am calculating energy and gradients externally.
Also, I will need the updated geometry after an optimization step but optimize seems to return only energy and wavefunction as outputs.

You sound confused about what the optimize function does. It does not take in a gradient and then return the next step. It takes in a function to compute gradients at a given geometry and runs the entire optimization from there. All the setting of geometries is handled internally. That’s why your current optimize call is raising an error. You need to send in a function to compute gradients, not the particular gradient returned by your grad_fun.get_grad call.

Thank you for the clarification. I misunderstood what optimize function was doing. The gradient calculations I am doing need not be Psi4 dependent, that’s why, I just need an optimizer, that is, which takes gradient as an input and returns the next step. Is there any other function in Psi4 which does what I need?

If that’s what you want, then your best option is to look at the GitHub repository that hokru posted earlier.

It seems like the optimizer in that GitHub repository can make use of different programs (psi4, pyscf, etc.) to compute gradients and energies. However, as I mentioned in my previous comments, I just need a function that takes in the gradients and provides an updated geometry of the molecule. Is there any function which does that in optking? It would be helpful if you could point me in the right direction.

Hi,
I apologize for how user-unfriendly this is. An interface for this is currently in planning.

I had to modify optkings __init__.py to get access to a method you would need. You can grab the version from my fork https://github.com/AlexHeide/optking.

This should do what you need using psi4 for easy molecule creation and easy creation of optkings molecule. Lemme know if there’s something else you need optking to do (or not to do).

import psi4
import optking
from optking import molsys, stepAlgorithms, hessian, displace, intcosMisc, history
from optking import optparams as op

h2o = psi4.geometry("""
        O
        H 1 1.0
        H 1 1.0 2 106.0
    """)

E = psi4.energy('mp2/cc-pvdz')
gradient = psi4.gradient('mp2/cc-pvdz').np.flatten()

molsys_obj = molsys.Molsys.from_psi4_molecule(h2o)
optking.make_internal_coords(molsys_obj)

# Same guess types as c++ optking
# intcosMisc has a convert_hessian_to_internals method if needed
H = hessian.guess(molsys_obj, guessType=op.Params.intrafrag_hess)
fq = intcosMisc.q_forces(molsys_obj.intcos, molsys_obj.geom, gradient)

# very user unfriendly but necessary sorry history object initialization
history.oHistory.append(molsys_obj.geom, E, fq)

print(molsys_obj.show_geom())

# None is energy, not needed for any algorithms except for the linesearch
step = stepAlgorithms.take_step(molsys_obj, None, fq, H, stepType='RFO')
# if you'd like to have optking do hessian updating

history.oHistory.hessian_update(H, molsys_obj)
print(molsys_obj.show_geom()) 

Please let me know if you run into any issues, and again I apologize that we don’t have a single method to do this for you.

I really appreciate you taking the time to add this feature in optking. I tried this script on a test case and it did provide me an updated geometry when I gave gradient and energy as the inputs. I didn’t run into any issues, but I should probably run more test cases to say that for sure.

I am glad to hear that an interface for external optimization is being planned for Psi4. I have used optimizers of Gaussian and Orca, externally. As a user, I would like the program to make the decision of when to stop the optimization procedure. The user should only be responsible for providing gradients and energy at each step of the optimization. For now, I think I can write a loop to perform the optimization and stop the procedure based on gradient norm or maximum gradient.

Thank you once again for your help.

Happy to help. An addition that would likely be helpful for your loop.

from optking import optwrapper, convCheck

optwrapper.initialize_options({"g_convergence": "foo_bar"})

If you would like to use optking’s convergence machinery you can call this method. Works as seen in the Psi4 optking documentation. q_pivot is for IRCs.

def convCheck.conv_check(iterNum, oMolsys, dq, f, energies, q_pivot=None):