Issues with IRC calculations following a frequency calculations

Hi!

I am running the following code to test the optimization of a transition state followed by frequency calculations and IRC test:

molecule ts {
0 1
 O  9.978538 -7.065826 -5.560229
 O  9.253519 -7.435897 -4.707331
 H  11.121586 -7.116893 -5.275842
 H  9.800901 -7.736312 -3.707248
} 

set {
	basis Def2-SVP
	reference rhf 	                # To ensure that the scf guess is right
	g_convergence gau               # Sets the SCF convergence as scf=tight in G16
	dft_spherical_points 302        # The two following options makes the grid as 'finegrid' in Gaussian
	dft_radial_points 75
	dft_pruning_scheme robust       # Generally safe and will speed things up
	opt_type ts   	                # Options min, ts, irc
	geom_maxiter 100
}

# First we optimize the TS
e, wfn = optimize('M06-2X-D3ZERO2B', return_wfn=True)
wfn.to_file('OPT')
ts.save_xyz_file('HOOH_ts_opt.xyz',1)

# Clone the final geometries for the IRC calculations
irc_forwards  = ts.clone()
irc_backwards = ts.clone()

# Now we get the frequencies
e, wfn_freq = frequency('M06-2X-D3ZERO2B', ref_gradient=wfn.gradient(), return_wfn=True)   # Save the frequency caluclation in the wavefunction for printing

# Now we move to the IRC caluclations
set {
	basis Def2-SVP
	reference rhf 	                # To ensure that the scf guess is right
	g_convergence gau               # Sets the SCF convergence as scf=tight in G16
	dft_spherical_points 302        # The two following options makes the grid as 'finegrid' in Gaussian
	dft_radial_points 75
	dft_pruning_scheme robust       # Generally safe and will speed things up
	opt_type irc   	                # Options min, ts, irc
	geom_maxiter 100
}

# IRC forward ####################################
set {
	irc_step_size 0.25
	irc_points 20
	irc_direction forward
	cart_hess_read = True
	full_hess_every = -1
}
e, irc_forw = optimize('M06-2X-D3ZERO2B', molecule=irc_forwards, ref_gradient=wfn_freq.gradient())
irc_forwards.save_xyz_file('ircf.xyz',1)

# IRC backward ####################################
set {
	irc_step_size 0.25
	irc_points 20
	irc_direction backward
	cart_hess_read = True
	full_hess_every = -1
}
e, irc_back = optimize('M06-2X-D3ZERO2B', molecule=irc_backwards, ref_gradient=wfn_freq.gradient())
irc_backwards.save_xyz_file('ircr.xyz',1)

Everything goes well for the optimization of the TS and the calculations of the frequencies.
My issues are two:

  1. If I add the “wfn_freq.to_file(‘wfn_freq’)” command after the frequency calculation (not shown in the code above), I get the same error I reported in a previous post here.

  2. If I skip that step and try to proceed directly with the IRC calculation as shown in the code above, I get the following error:

Printing out the relevant lines from the Psithon --> Python processed input file:
    core.set_global_option("IRC_STEP_SIZE", 0.25)
    core.set_global_option("IRC_POINTS", 20)
    core.set_global_option("IRC_DIRECTION", "forward")
    core.set_global_option("CART_HESS_READ", "True")
    core.set_global_option("FULL_HESS_EVERY", -1)
--> e, irc_forw = optimize('M06-2X-D3ZERO2B', molecule=irc_forwards, ref_gradient=wfn_freq.gradient())
    irc_forwards.save_xyz_file('ircf.xyz',1)
    core.set_global_option("IRC_STEP_SIZE", 0.25)
    core.set_global_option("IRC_POINTS", 20)
    core.set_global_option("IRC_DIRECTION", "backward")
    core.set_global_option("CART_HESS_READ", "True")

!-----------------------------------------------------------------------!
!                                                                       !
!  [Errno 2] No such file or directory: 'HOOH_freq_irc.ts.1706054.hess' !
!                                                                       !
!-----------------------------------------------------------------------!

Hi!

For this, all you need to do is add set hessian_write True before the frequency calculation. It also looks like you’ll want to add a return_wfn=True to your irc optimizations.

1 Like

Thanks!
I placed that command just before the frequency calculation but unfortunately I still get the other error I flagged in the other post I linked here:

Printing out the relevant lines from the Psithon --> Python processed input file:
    ts.save_xyz_file('HOOH_ts_opt.xyz',1)
    irc_forwards  = ts.clone()
    irc_backwards = ts.clone()
    core.set_global_option("HESSIAN_WRITE", "True")
    e, wfn_freq = frequency('M06-2X-D3ZERO2B', ref_gradient=wfn.gradient(), return_wfn=True)  
--> wfn_freq.to_file('FREQ')
    core.set_global_option("BASIS", "Def2-SVP")
    core.set_global_option("REFERENCE", "rhf")
    core.set_global_option("G_CONVERGENCE", "gau")
    core.set_global_option("DFT_SPHERICAL_POINTS", 302)
    core.set_global_option("DFT_RADIAL_POINTS", 75)

!----------------------------------------------------------------------------------!
!                                                                                  !
! Fatal Error: Wavefunction::Ca: Unable to obtain MO coefficients.                 !
! Error occurred in file: /scratch/psilocaluser/conda-                             !
!     builds/psi4-multiout_1670980379431/work/psi4/src/psi4/libmints/wavefunction. !
!     cc on line: 804                                                              !
! The most recent 5 function calls were:                                           !
! psi::Wavefunction::Ca() const                                                    !
!                                                                                  !
!----------------------------------------------------------------------------------!

Update:

My bad, in the previous error I still had the line to save the wfn from the frequency caluclation:

e, wfn_freq = frequency('M06-2X-D3ZERO2B', ref_gradient=wfn.gradient(), return_wfn=True)
wfn_freq.to_file('FREQ')

but if I remove it then the calculation proceeds perfectly and I see that the Hessian is saved as a .hess.
I guess that if I want to perform the two calculations, frequency and IRC, in different times I can recall the Hessian file in the IRC input, right? (to not re-calculate frequencies from scratch again).

Correct. If you want to run a new calculation with a .hess you previously generated you’ll just need to rename it with something like.

import os
os.system(f"cp HOOH_freq_irc.ts.1706054.hess HOOH_freq_irc.ts.{os.getpid()}.hess")

This just properly associates it with the current psi4 process.

That great! But I have some issues on how to recall the .hess file for the IRC calculation (not sure what the keyword is and I cannot find it on the doc’s).
I have all the options for the IRC caluclations set but how do I “link” the previous .hess file to the IRC calculation; i.e., should I still use the ref_gradient function? If you could fill the gaps in the block below it would be great!

import os
os.system(f"cp HOOH_freq_irc.ts.1706054.hess HOOH_freq_irc.ts.{os.getpid()}.hess")

# Probably need to assign the new .hess with a variable name I guess and actually load the .hess file as well

# Here I don't know how to tell the IRC calculation to use the previously loaded .hess file to start the IRC process.
e, irc_forw = optimize('M06-2X-D3ZERO2B', molecule=irc_forwards, return_history=True, ref_gradient=XXXXXXX)


Hi!

All you should need to do is copy the .hess file so that it has the correct pid and run the irc with
cart_hess_read True. What you had in your original input file should have worked:

set {
	basis Def2-SVP
	reference rhf 	                # To ensure that the scf guess is right
	g_convergence gau               # Sets the SCF convergence as scf=tight in G16
	dft_spherical_points 302        # The two following options makes the grid as 'finegrid' in Gaussian
	dft_radial_points 75
	dft_pruning_scheme robust       # Generally safe and will speed things up
	opt_type irc   	                # Options min, ts, irc
	geom_maxiter 100
}

# IRC forward ####################################
set {
	irc_step_size 0.25
	irc_points 20
	irc_direction forward
	cart_hess_read = True
	full_hess_every = -1
}
e, irc_forw = optimize('M06-2X-D3ZERO2B', molecule=irc_forwards, ref_gradient=wfn_freq.gradient())
irc_forwards.save_xyz_file('ircf.xyz',1)

You don’t need to pass a reference gradient. The ref_gradient argument is only utilized afaik in frequency calculations to determine if the finite difference calculation is being performed at a stationary point.

Fantastic! Thanks for the clarification - I better understand the logic behind it now.
I tried it in a couple of systems and it seems to work fine.
I have also tried it on a transition state for a bisradical species and I got this error:

Traceback (most recent call last):
  File "/home/ivan/mambaforge/envs/psi4_env/bin/psi4", line 338, in <module>
    exec(content)
  File "<string>", line 97, in <module>
  File "/home/ivan/mambaforge/envs/psi4_env/lib//python3.9/site-packages/psi4/driver/driver.py", line 1272, in optimize
    opt_object.take_step()
  File "/home/ivan/mambaforge/envs/psi4_env/lib//python3.9/site-packages/optking/opt_helper.py", line 174, in take_step
    self.dq, self.step_str = self.opt_manager.take_step(self.fq, self._Hq, self.E, return_str=True)
  File "/home/ivan/mambaforge/envs/psi4_env/lib//python3.9/site-packages/optking/optimize.py", line 232, in take_step
    achieved_dq, returned_str = self.opt_method.take_step(fq, H, energy, return_str=True)
  File "/home/ivan/mambaforge/envs/psi4_env/lib//python3.9/site-packages/optking/IRCfollowing.py", line 97, in take_step
    dq = self.dq_irc(fq, H)
  File "/home/ivan/mambaforge/envs/psi4_env/lib//python3.9/site-packages/optking/IRCfollowing.py", line 212, in dq_irc
    G_prime_root_inv = symm_mat_root(G_prime_inv)
  File "/home/ivan/mambaforge/envs/psi4_env/lib//python3.9/site-packages/optking/linearAlgebra.py", line 171, in symm_mat_root
    rootMatrix[i][i] = sqrt(evals[i])

ValueError: math domain error


Printing out the relevant lines from the Psithon --> Python processed input file:
    core.set_global_option("IRC_STEP_SIZE", 0.25)
    core.set_global_option("IRC_POINTS", 20)
    core.set_global_option("IRC_DIRECTION", "forward")
    core.set_global_option("CART_HESS_READ", "True")
    core.set_global_option("FULL_HESS_EVERY", -1)
--> e, irc_forw = optimize('M06-2X-D3ZERO2B', molecule=irc_forwards, return_wfn=True)
    irc_forwards.save_xyz_file('TS_radH_abstraction_primary_00_ircf.xyz',1)
    core.set_global_option("IRC_STEP_SIZE", 0.25)
    core.set_global_option("IRC_POINTS", 20)
    core.set_global_option("IRC_DIRECTION", "backward")
    core.set_global_option("CART_HESS_READ", "True")

!--------------------------------!
!                                !
!  math domain error             !
!                                !
!--------------------------------!

    Psi4 stopped on: Thursday, 06 April 2023 05:43PM
    Psi4 wall time for execution: 0:06:13.14

*** Psi4 encountered an error. Buy a developer more coffee!
*** Resources and help at github.com/psi4/psi4.

Do you think that this is something related to the system/PSI4 or something else?

Hi Fransu.
For a question like this, can you upload the input and output files? The G matrix should be positive definite. I wonder if something may be wrong or nonphysical with the geometry for such an error to occur.