IRC Step Back-Transformation Not Converging

I have been trying to run an IRC using the latest version of Psi4 (specifically commit 02a1496), but the IRC fails with the following message:

Read in cartesian Hessian and transformed it.

IRC_DATA is empty, so we are at the transition state.
Stepping in forward direction from TS.

Back-transformation to cartesian coordinates...
Could not converge backtransformation.
Using first guess instead.
The INTCO_EXCEPTion handler:
Could not take constrained step in an IRC computation.
Dynamic level is 0.
exc.g_really_quit() is 0.

**** Optimization has failed! (in 1 steps) ****

Is there anything I can do to make this converge?

A few things that may help:

  • Because of unrelated technical difficulties, it is more convenient for me to read in a Hessian produced by a separate job than it is to compute an initial Hessian here. I have tried setting write_hess_every to 0 instead. I get the same problem. If you want, I can supply you with the Hessian I’ve been using.
  • I have the same problem whether I run the IRC in the backwards or the forwards direction.
  • While I have not had difficulties converging this system in CFOUR previously, the expected reactant is a hydrogen-bonded complex.

Input file below:

memory 28 Gb

molecule CI {
symmetry c1
O
O 1 r1
C 2 r2 1 a1
H 3 r3 2 a2 1 d1
H 3 r8 2 a7 1 d2
H 1 r5 2 a4 3 d4
N 3 r4 6 a3 1 d3
H 7 r6 6 a5 1 d5
H 7 r7 6 a6 8 d6

r1   =        1.480966039149124
r2   =        1.271341026926654
a1   =      108.047828274256858
r3   =        1.092332520584757
a2   =      116.065939142851818
d1   =      171.596874434455515
r8   =        1.088922891150393
a7   =      119.480789128455314
d2   =      -29.742352341603691
r5   =        1.878515775552397
a4   =       81.905285481429161
d4   =      -53.152262596000298
r4   =        2.089375187124503
a3   =       28.068484657311583
d3   =     -176.996973742810468
r6   =        1.018034441266780
a5   =      109.865405959970531
d5   =      129.409401731653560
r7   =        1.018450111243980
a6   =      109.181643309544015
d6   =      119.598212102759803
}

set {
basis ano0
opt_type irc
g_convergence gau_verytight
irc_direction forward
cart_hess_read true
writer_file_label CI
geom_maxiter 150
print_trajectory_xyz_file true
irc_step_size 0.05
intrafrag_step_limit 0.05
intrafrag_step_limit_max 0.05
}

optimize('c4-ccsd(t)')

I’m testing this at a lower level of theory on my system, but maybe @Rollin_King has some insights off the bat?

I’ve done some additional testing at the SCF level, and I’ve gotten some strange results. I can get convergence if I use the ensure_bt_convergence keyword, but the step size never changes. The first time you attempt convergence with the keyword enabled, fix_bend_axes is not called. If you do not include the keyword but comment out the fix_bend_axes call on line 149 here, you get convergence.

Including the keyword also gets CCSD(T) convergence.

The keyword “ensure_bt_convergence” reduces the step size until one can be taken in cartesians. However, the first event in our IRC algorithm is a step in the direction of a Hessian eigenvector that is 1/2 the length of the requested IRC step size. So this keyword is probably not helping you get what you need.

On the other hand, if commenting out ‘fix_bend_axes’ allows convergence then I believe you are OK. This keyword fixes the meaning of quasi-linear bends. In some cases, it is needed to get a successful back-transformation but in others I’ve recently found it prevents a precisely chosen step. The latter is not generally critical for stationary point searches.

It’s a good idea to test with SCF if you can. Also, with printing turned up, the Hessian after its transformation to internal coordinates will be printed in the output. You could visually inspect that to see if it is reasonable. That is, do you see larger diagonal elements appropriate for stretches and bends.

At the SCF level, my Hessian looks reasonable, and I’m getting a reasonable trajectory. It looks like this was just one of those cases where “fix_bend_axes” shouldn’t be used, thanks!

We should probably modify the code to try to accomplish the step without
fixing the axes, and then if that fails (might be tricky to define “fail”
in general) retry with the axes fixed.

Separate issue – for reading in a Hessian that wasn’t computed within the same job, I believe you are going to run into problems. Specifically, at some point the libmints function get_writer_file_prefix began including the process ID – it wasn’t a problem before that. Since there is no way to know the PID in advance, this prevents the user from reading in a Hessian that was computed in a separate job. I think it might be worth avoiding the use of get_writer_file_prefix within optking and allowing the user to set a file prefix for optking, or perhaps just for reading the Hessian. Is that right @Rollin_King – what do you think?

A temporary work-around would be to rename the Hessian file within your IRC job input using the current PID, but that seems seriously sub-optimal as a long-term solution.

Looks like this change was made back in January for MPI compatibility.

Not sure what the considerations are there, as far as OPTKING is concerned – I don’t think OPTKING uses MPI.

I can get around this because my cluster manager modifies Psi4 so it can read in the Hessian without a PID. This, of course, does nothing for the general case.

Some users have reported problems with supplying files due to a lack of a PID here. One of these files is the internal coordinate file for an optimization.

Yeah, we need a hack for this. There was, once upon a time, a PID override
that could be provided by the user, but I believe it has been obsoleted.
As long as we are using PIDs in filenames, then I think this functionality
should be brought back. There will always be at least ‘expert’ uses for it.