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.core.clean_options()
psi4.set_memory("4GB")
psi4.set_output_file("cyclobutene_opening_IRC_10.log", append=False)
psi4.core.set_num_threads(4)
### 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
psi4.set_options({
"geom_maxiter":500,
"full_hess_every":1, # Tried -1 and 0
"opt_type":"irc",
"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
})
print("start")
E, history = psi4.optimize('b3pw91/6-31+G(d)',molecule=mol, return_history=True)
print("done")
The calculation seems to complete two steps and then fails with the following output generated…
start
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/driver.py:1272, 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/opt_helper.py:174, 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/optimize.py:232, 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/IRCfollowing.py:97, 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 logger.info("IRC Constrained step calculation finished.")
File ~/psi4conda/lib/python3.11/site-packages/optking/IRCfollowing.py:212, 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/linearAlgebra.py:171, 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 = np.dot(evects, np.dot(rootMatrix, 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?