Geometry Optimisation with arbitrary wavefunction

feature-request
1

Dear Psi4 team,

I’ve been using your program (python interface) to build an SCF method with corrections for dispersion added to the core Hamiltonian, with good results. We’d like to move to performing geometry optimisations, but this is notably limited to the built in methods.

I note that the optimise function works by repeatedly calling gradients and OPTKING. Assuming I can generate the gradients for the additional terms independently (analytically/numerically/whatever), is there a straightforward way to add them to the optimisation step?

Cheers,
Jed

2

This is something I’ve needed for my own research! I’m going to warn you that this part of the interface is going to be changing soon as part of a larger change to the driver logic and the optimizer.

Here’s what you can do for now:
Call psi4.optimize(your-function-here, dertype=1, molecule=your-molecule-here). The signature Psi expects of your function is shown here. gradient is Psi’s own gradient method, which is passed as the first argument for historical reasons. Your function must return a Psi4 Matrix object as the first argument and a second argument. Although Psi expects this to be a wavefunction, nothing bad will happen if you have something else as your second argument, as long as you don’t pass return_wfn=True to your psi4.optimize call.

In your case, your gradient function will probably be a regular Psi4 gradient call, a call to your dispersion function, and then returning the sum of the two.

3

Thanks Jonathan, appreciate the quick reply. I’ll give it a try and see how we go.

Great minds on the research :wink:

Cheers,
Jed

4

This is a hack that will not work at some unspecified time in the future. But I thought it worth documenting that if you don’t need the optimize() bells and whistles (finite differences, occasional Hessian compuations, etc.) then you can hardwire an optimization like this. This is a toy example of a harmonic force on a couple atoms, but you’ll get the idea.

molecule n2 {
  0 1
  N  0.0  0.0 -1.0
  N  0.0  0.0  1.0
}
set basis cc-pvdz

mol = psi4.get_active_molecule()
mol.update_geometry()
core.set_legacy_molecule(mol)

k = 0.1

for step in range(20):
    xyz = mol.geometry().np

    grad = np.array( [
      [0.0, 0.0, 0.5*k*(xyz[0,2]-1.3)],
      [0.0, 0.0, 0.5*k*(xyz[1,2]+1.3)] ] )
    core.set_legacy_gradient(Matrix.from_array(grad))

    E  = k*(xyz[0,2]-1.3)*(xyz[0,2]-1.3) + k*(xyz[1,2]+1.3)*(xyz[1,2]+1.3)
    core.set_variable('CURRENT ENERGY', E)

    rval = core.optking()
    if rval != 0:
        break
5

Better to switch over to the new python optimizer, which is not hard to do, but will be easy in the near future.