IRC optimization fails with "Math Domain Error"

I have been building Jupyter notebooks using Psi4 for my Physical organic Chemistry class. Everything was going well until I tried to build a reaction coordinate using IRC optimization.

I am using an IRC optimization of a transition state for thermal con-rotatory ring opening of cyclobutene. I was able to optimize the ts and confirm the imaginary vibration. Then I set up the following Python code to perform an IRC calculation…

# use psi4conda environment
import psi4
import numpy as np

psi4.set_output_file("cyclobutene_opening_IRC_10.log", append=False)

### optimized transition state structure
data = """
       0 1
 C    1.040389109196    0.254964976384   -0.647674881834
 C    0.689591079928    0.000000000000    0.708001346062
 C   -0.689591079928    0.000000000000    0.708001346062
 C   -1.040389109196   -0.254964976384   -0.647674881834
 H    1.923694934595   -0.185448301950   -1.124148093391
 H    0.711123863746    1.180002035639   -1.108938410294
 H    1.381683913454   -0.238567611354    1.514789629143
 H   -1.381683913454    0.238567611354    1.514789629143
 H   -1.923694934595    0.185448301950   -1.124148093391
 H   -0.711123863746   -1.180002035639   -1.108938410294

       units angstrom
mol = psi4.geometry(data)             # Create Molecule object from data string

     "full_hess_every":1,                # Tried -1 and 0
     "irc_step_size":0.15,               # Also tried 0.20 and 0.25
     "ensure_bt_convergence":True,       # Tried False
     "irc_direction": "forward",         # Tried 'backward'
     "reference": "rhf",                 # Also tried "uhf"
     "g_convergence": "GAU_VERYTIGHT",   # Also tried GAU_TIGHT and QCHEM
     "IRC_POINTS": 10                     # Tried 20 and 25
E, history = psi4.optimize('b3pw91/6-31+G(d)',molecule=mol, return_history=True)

The calculation seems to complete two steps and then fails with the following output generated…

Warning: thermodynamics relations excluded imaginary frequencies: ['732.1013i']
Warning: used thermodynamics relations inappropriate for low-frequency modes: ['472.1911']
Warning: thermodynamics relations excluded imaginary frequencies: ['722.9255i']
Warning: used thermodynamics relations inappropriate for low-frequency modes: ['474.7620']
ValueError                                Traceback (most recent call last)
Cell In[2], line 56
     44 psi4.set_options({
     45      "geom_maxiter":500,
     46      "full_hess_every":1,                # Tried -1 and 0
     53      "IRC_POINTS": 5                     # Tried 20 and 25
     54      })
     55 print("start")
---> 56 E, history = psi4.optimize('b3pw91/6-31+G(d)',molecule=mol, return_history=True)
     57 print("done")

File ~/psi4conda/lib/python3.11/site-packages/psi4/driver/, in optimize(name, **kwargs)
   1270 opt_object.compute()  # process E, gX, H
   1271 try:
-> 1272     opt_object.take_step()
   1273 except optking.exceptions.AlgError:
   1274     # Optking encountered an algorithm error and reset.
   1275     if not opt_object.HX:

File ~/psi4conda/lib/python3.11/site-packages/optking/, in Helper.take_step(self)
    172 self.opt_manager.error = None
    173 try:
--> 174     self.dq, self.step_str = self.opt_manager.take_step(self.fq, self._Hq, self.E, return_str=True)
    175 except AlgError as e:
    176     self.opt_manager.alg_error_handler(e)

File ~/psi4conda/lib/python3.11/site-packages/optking/, in OptimizationManager.take_step(self, fq, H, energy, return_str, **kwargs)
    229 self.current_requirements = self.update_requirements()
    231 if not self.params.linesearch:
--> 232     achieved_dq, returned_str = self.opt_method.take_step(fq, H, energy, return_str=True)
    233 else:
    234     if self.direction is None:

File ~/psi4conda/lib/python3.11/site-packages/optking/, in IntrinsicReactionCoordinate.take_step(self, fq, H, energy, return_str, **kwargs)
     94 else:
     96     self.history.append(self.molsys.geom, energy, fq, self.molsys.gradient_to_cartesians(-1 * fq))
---> 97     dq = self.dq_irc(fq, H)
     98     dq, dx, return_str = displace_molsys(self.molsys, dq, fq, ensure_convergence=True, return_str=True)
     99"IRC Constrained step calculation finished.")

File ~/psi4conda/lib/python3.11/site-packages/optking/, in IntrinsicReactionCoordinate.dq_irc(self, f_q, H_q)
    210 G_prime_root = symm_mat_root(G_prime)
    211 G_prime_inv = symm_mat_inv(G_prime, redundant=True)
--> 212 G_prime_root_inv = symm_mat_root(G_prime_inv)
    214 logger.debug("G prime root matrix: \n" + print_mat_string(G_prime_root))
    216 deltaQM = 0

File ~/psi4conda/lib/python3.11/site-packages/optking/, in symm_mat_root(A, Inverse)
    168         evals[i] = 1 / evals[i]
    170 for i in range(0, len(evals)):
--> 171     rootMatrix[i][i] = sqrt(evals[i])
    173 A =,, evects.T))
    175 return A

ValueError: math domain error

I have taken the optimized transition state geometry that is in the input file above and used it in an IRC calculation in GAMESS using the same level of theory and it worked perfectly. So I don’t believe it is an issue with the geometry or method.

I have uploaded the output file. I have run the code in a Jupyter notebook and as a stand-alone python program with the same results.

cyclobutene_opening_IRC_10.log.txt (1.7 MB)

I note the mention of a "Math Domain Error’ in a recent thread…

… but no resolution was reached for that error in that discussion.

I have tried many different options and can’t get past this error. Any suggestions?

UPDATE: I tried the data from the two IRC test files on the Psi4 GitHub site: psi4/samples/opt-irc-1/input.dat at master · psi4/psi4 · GitHub
psi4/samples/opt-irc-2/input.dat at master · psi4/psi4 · GitHub

One is for the IRC starting from the TS for rotation in HOOH and the other is for the migration of H in HCN.

Both worked. I then tried the code in the question above again and git the same “Math Domain Error” failure.

Below is the python version of the HCN test that worked. I used the same setup as in the ring-opening IRC that gives the error.

# use psi4conda environment
import psi4


# The Z-matrix as a text string ### From previous optimization
data = """
       0 1
  H    -0.5958806528   0.9889214459   0.0000000000
  C    -0.5958806528  -0.1660941336   0.0000000000
  N     0.5535292657   0.0711607905   0.0000000000
       units angstrom

# Create the Molecule object
mol = psi4.geometry(data)             # Create Molecule object from data string

output_file = "HCN_test.log"
psi4.set_output_file(output_file, append=False, loglevel=20, print_header=True, inherit_loglevel=True, execute=True)

irc_direction = "forward"

     "reference": "rhf",
     "IRC_POINTS": 10

E, history = psi4.optimize('b3pw91/6-31+G(d)',molecule=mol, return_history=True)


So the mystery is that the TS for cyclobutene ring opening is able to produce a successful IRC in GAMESS but in Psi4 it fails with the “Math Domain Error” described above.

I will keep poking away at this. Any advice would be appreciated.