Frequency with PCM takes forever for finite difference calculation

Hi,
I am very new to PSI4 and also I am not really a pyhsical chemist and trying to reproduce some benchmark tests from a paper that my team is interested in.

I am testing both PSI4 and Jaguar but somehow, PSI4 takes almost 24 hours for one single frequency calculation while Jaguar takes only 30 minutes with the same basis set and resources. I don’t believe PSI4 is supposed to be this slow.
Especially, it gets stuck with this message “140 displacements needed.”
Would this related to the symmetry that I set in the geometry? PCM forces user to specify the symmetry and I was not sure what to put…

Also, when I look at the CPU usage, it takes the full CPU without PCM while it takes only one thread with PCM. Is it because PCM in PSI4 is not opimized/parallelized?

Can anyone please see what is wrong with my script?

import os
import psi4
psi4.core.set_output_file(‘output.dat’, False)
psi4.core.set_num_threads(12)
psi4.set_memory(‘16 GB’)
print(‘You are using {} threads.’.format(psi4.core.get_num_threads() ))

mol = psi4.geometry("""
symmetry c1
C 4.4314 1.0413 1.4383
S 4.5445 -0.3623 0.2710
C 2.7974 -0.9102 0.2525
C 1.8404 0.0654 -0.4391
C 0.4350 -0.5302 -0.5327
O 0.2521 -1.7226 -0.7391
N -0.5779 0.3888 -0.3818
C -1.9693 0.2024 -0.4402
C -2.5810 -1.0373 -0.6829
N -2.6479 1.3397 -0.2362
C -3.9712 -1.0679 -0.7092
C -3.9829 1.2789 -0.2672
C -4.6979 0.1055 -0.4984
H 5.4516 1.3928 1.6148
H 4.0005 0.7206 2.3927
H 3.8448 1.8703 1.0307
H 2.4718 -1.1061 1.2808
H 2.7811 -1.8664 -0.2767
H 1.8256 1.0310 0.0795
H 2.1933 0.2574 -1.4604
H -0.3320 1.3550 -0.2038
H -1.9765 -1.9193 -0.8393
H -4.4829 -2.0086 -0.8939
H -4.5017 2.2207 -0.0973
H -5.7831 0.1127 -0.5124
units Angstrom
“”")

print(“Optimizing geometry…”)
psi4.optimize(‘b3lyp-d3/6-31G(d)’, molecule=mol)

psi4.set_options({
‘basis’:‘6-311+G(d,p)’,
‘scf_type’:‘df’,
‘pcm’: True,
‘pcm_scf_type’: ‘total’,
})

psi4.pcm_helper("""
Units = Angstrom
Medium {
SolverType = IEFPCM
Solvent = Water
}
Cavity {
RadiiSet = UFF
Type = GePol
Scaling = False
Area = 0.3
Mode = Implicit
}
“”")

print(‘Calculating single point energy…’)
E, wfn = psi4.frequency(‘b3lyp-d3/6-311+G(d,p)’, return_wfn=True)

Unfortunately only rudimentary PCM feature are available at this time.
At least part of the PCM part is not threaded, and more importantly no analytical gradients are available.
Double numerical derivatives (numerical hessian with numerical gradient) are often unstable.
PCM cannot make use of symmetry, you need to specify c1.

There are plans to get gradients for implicit solvation models, likely with a different PCM library, but that is not expected in the near future.

Short term better integral screening will lead to faster PCM calculations, but frequencies will remain very tricky.