Frequency calculations result in `<built-in method X of PyCapsule object at Y> returned NULL without setting an error`

Hi,

I’m experiencing strange behavior when conducting frequency calculations on two different geometries. I’m using the PBE0method with the cc-pvdz basis set, and when calling psi4.frequency(), I get the following errors for the optimized H2O and HOCl geometries:

Traceback for HOCl:

~/anaconda3/envs/psi4/lib/python3.6/site-packages/psi4/driver/driver.py in frequency(name, **kwargs)
   1473     # Compute the hessian
-> 1474     H, wfn = hessian(name, return_wfn=True, molecule=molecule, **kwargs)
   1475 

~/anaconda3/envs/psi4/lib/python3.6/site-packages/psi4/driver/driver.py in hessian(name, **kwargs)
   1214     elif gradient_type == "cbs_gufunc":
-> 1215         return driver_cbs._cbs_gufunc(hessian, name.lower(), **kwargs, ptype="hessian")
   1216     elif gradient_type == "cbs_wrapper":

~/anaconda3/envs/psi4/lib/python3.6/site-packages/psi4/driver/driver_cbs.py in _cbs_gufunc(func, total_method_name, **kwargs)
   1959         core.set_global_option('BASIS', basis)
-> 1960         ptype_value, wfn = func(method_name, return_wfn=True, molecule=molecule, **kwargs)
   1961         core.clean()

~/anaconda3/envs/psi4/lib/python3.6/site-packages/psi4/driver/driver.py in hessian(name, **kwargs)
   1266     else:
-> 1267         G0 = gradient(lowername, molecule=molecule, **kwargs)
   1268     translations_projection_sound, rotations_projection_sound = _energy_is_invariant(G0)

~/anaconda3/envs/psi4/lib/python3.6/site-packages/psi4/driver/driver.py in gradient(name, **kwargs)
    698         # Perform the gradient calculation
--> 699         wfn = procedures['gradient'][lowername](lowername, molecule=molecule, **kwargs)
    700 

~/anaconda3/envs/psi4/lib/python3.6/site-packages/psi4/driver/procrouting/proc.py in run_scf_gradient(name, **kwargs)
   2139     if ref_wfn is None:
-> 2140         ref_wfn = run_scf(name, **kwargs)
   2141 

~/anaconda3/envs/psi4/lib/python3.6/site-packages/psi4/driver/procrouting/proc.py in run_scf(name, **kwargs)
   2044 
-> 2045     scf_wfn = scf_helper(name, post_scf=False, **kwargs)
   2046     returnvalue = scf_wfn.energy()

~/anaconda3/envs/psi4/lib/python3.6/site-packages/psi4/driver/procrouting/proc.py in scf_helper(name, post_scf, **kwargs)
   1362 
-> 1363     e_scf = scf_wfn.compute_energy()
   1364     for obj in [core, scf_wfn]:

~/anaconda3/envs/psi4/lib/python3.6/site-packages/psi4/driver/procrouting/scf_proc/scf_iterator.py in scf_compute_energy(self)
     82     else:
---> 83         self.initialize()
     84 

~/anaconda3/envs/psi4/lib/python3.6/site-packages/psi4/driver/procrouting/scf_proc/scf_iterator.py in scf_initialize(self)
    196 
--> 197         mints.one_electron_integrals()
    198         self.initialize_jk(self.memory_jk_, jk=jk)

SystemError: <built-in method one_electron_integrals of PyCapsule object at 0x7fe08c357750> returned NULL without setting an error

Traceback for H2O:

~/anaconda3/envs/psi4/lib/python3.6/site-packages/psi4/driver/driver.py in frequency(name, **kwargs)
   1473     # Compute the hessian
-> 1474     H, wfn = hessian(name, return_wfn=True, molecule=molecule, **kwargs)
   1475 

~/anaconda3/envs/psi4/lib/python3.6/site-packages/psi4/driver/driver.py in hessian(name, **kwargs)
   1214     elif gradient_type == "cbs_gufunc":
-> 1215         return driver_cbs._cbs_gufunc(hessian, name.lower(), **kwargs, ptype="hessian")
   1216     elif gradient_type == "cbs_wrapper":

~/anaconda3/envs/psi4/lib/python3.6/site-packages/psi4/driver/driver_cbs.py in _cbs_gufunc(func, total_method_name, **kwargs)
   1959         core.set_global_option('BASIS', basis)
-> 1960         ptype_value, wfn = func(method_name, return_wfn=True, molecule=molecule, **kwargs)
   1961         core.clean()

~/anaconda3/envs/psi4/lib/python3.6/site-packages/psi4/driver/driver.py in hessian(name, **kwargs)
   1294         # Obtain list of displacements
-> 1295         findif_meta_dict = driver_findif.hessian_from_gradients_geometries(molecule, irrep)
   1296 

~/anaconda3/envs/psi4/lib/python3.6/site-packages/psi4/driver/driver_findif.py in hessian_from_gradients_geometries(molecule, irrep)
    888     
--> 889     return _geom_generator(molecule, irrep, "2_1")
    890 

~/anaconda3/envs/psi4/lib/python3.6/site-packages/psi4/driver/driver_findif.py in _geom_generator(mol, freq_irrep_only, mode)
    341 
--> 342     data = _initialize_findif(mol, freq_irrep_only, mode, init_string, 1)
    343 

~/anaconda3/envs/psi4/lib/python3.6/site-packages/psi4/driver/driver_findif.py in _initialize_findif(mol, freq_irrep_only, mode, initialize_string, verbose)
    170     # Populate salc_indices_pi for all irreps.
--> 171     for i, salc in enumerate(salc_list):
    172         salc_indices_pi[salc.irrep_index()].append(i)

SystemError: <built-in method __next__ of PyCapsule object at 0x7fe030db2f60> returned NULL without setting an error

I am able to optimize the two geometries just fine, and the logs don’t show any errors. They just cut out at:

H2O, end of log file:

  Based on options and gradient (rms=3.16E-05), recommend projecting translations and projecting rotations.
hessian() will perform frequency computation by finite difference of analytic gradients.

         ----------------------------------------------------------
                                   FINDIF
                     R. A. King and Jonathon Misiewicz
         ---------------------------------------------------------

  Using finite-differences of gradients to determine vibrational frequencies and 
  normal modes. Resulting frequencies are only valid at stationary points.
    Generating geometries for use with 3-point formula.
    Displacement size will be 5.00e-03.
    Number of atoms is 3.
    Number of irreps is 2.
    Number of SALCs is 3.
    Translations projected? 1. Rotations projected? 1.

HOCl, end of log file:

  ==> Pre-Iterations <==

   -------------------------------------------------------
    Irrep   Nso     Nmo     Nalpha   Nbeta   Ndocc  Nsocc
   -------------------------------------------------------
     A'        27      27       0       0       0       0
     A"        10      10       0       0       0       0
   -------------------------------------------------------
    Total      37      37      13      13      13       0
   -------------------------------------------------------

  ==> Integral Setup <==

I experienced both of these errors on both the latest v1.3 build and the latest dev build. Any ideas on what is going wrong?

Thanks.

Can you give input files so we can reproduce these?

Whups, sorry I forgot to provide that. I’m using the Python API, so these were my calls:

psi4_mol = """O   -1.97217300000000e-06   6.74217598520000e-02   0.00000000000000e+00
H    7.53506313381000e-01  -5.34994165006000e-01   0.00000000000000e+00
H   -7.53475013567000e-01  -5.35038090449000e-01   0.00000000000000e+00
0 1"""
total_energy, wavefunction = psi4.frequency("{}/{}".format('PBE0', 'cc-pvdz'), molecule=psi4_mol, return_wfn=True)

psi4_mol = """CL  -5.63286174164000e-01   2.40406396720000e-02   0.00000000000000e+00
O    1.14265757897400e+00  -1.04794065702000e-01   0.00000000000000e+00
H    1.40972964713600e+00   8.29011504722000e-01   0.00000000000000e+00
0 1"""
total_energy, wavefunction = psi4.frequency("{}/{}".format('PBE0', 'cc-pvdz'), molecule=psi4_mol, return_wfn=True)

Can you double-check the version of Psi you’re using? I’m using 1.3.2 release and get an error message due to frequency wanting an actual Molecule object, not a string, as the molecule argument.

If you’re sure that runs on 1.3.2, I can double check if something is amiss on my version.

I am running the latest conda version of 1.3.2. Also, my bad, I transferred the code wrong! I’m actually using a set of libraries that the lab I’m working in built as an automated wrapper around psi4, so I’ve tried to transfer the code manually with printed geometry strings and log statements and the like. It should look like this:

psi4_mol = psi4.geometry("""O   -1.97217300000000e-06   6.74217598520000e-02   0.00000000000000e+00
H    7.53506313381000e-01  -5.34994165006000e-01   0.00000000000000e+00
H   -7.53475013567000e-01  -5.35038090449000e-01   0.00000000000000e+00
0 1""")
total_energy, wavefunction = psi4.frequency("{}/{}".format('PBE0', 'cc-pvdz'), molecule=psi4_mol, return_wfn=True)

psi4_mol = psi4.geometry("""CL  -5.63286174164000e-01   2.40406396720000e-02   0.00000000000000e+00
O    1.14265757897400e+00  -1.04794065702000e-01   0.00000000000000e+00
H    1.40972964713600e+00   8.29011504722000e-01   0.00000000000000e+00
0 1""")
total_energy, wavefunction = psi4.frequency("{}/{}".format('PBE0', 'cc-pvdz'), molecule=psi4_mol, return_wfn=True)

Although, strangely enough, when I execute that code, it seems to work fine. I’ll dig a little deeper into my setup to try to see what differences are present.

Okay, I think I’ve found the source of the error. It seems to appear only when I run the code from within a Jupyter Notebook. When I run both the code snippet and my wrapper library from the Python shell or from an executable .py file, it works fine, but I get the PyCapsule errors when running some functions (I’ve experienced this error before with energy calculations as well, but I’ve been unable to reproduce those since it works for some geometries) in Jupyter Noteboks.

Are you sure that Jupyter and the command line use the same version of Python and Psi?

I am certain; I ran the Jupyter notebook from the conda environment in which I installed psi4. I do not have psi4 installed in my base environment, either.

I’ve learned to be paranoid for errors involving different ways of calling Python. Would you mind humoring me?

Run import psi4; psi4.__file__ from both the command line python and Jupyter, as well as import sys; sys.executable. It should give the same result in both ways of calling Python.

Sorry about the delay, I was away from my computer over the weekend.

Here’s the output in the Jupyter notebook:

import psi4
print(psi4.__file__)

import sys
print(sys.executable)

/server-home1/rroy/anaconda3/envs/psi4/lib/python3.6/site-packages/psi4/__init__.py
/server-home1/rroy/anaconda3/envs/psi4/bin/python

and here’s the output in the Python shell:

>>> import psi4
>>> psi4.__file__
'/server-home1/rroy/anaconda3/envs/psi4/lib/python3.6/site-packages/psi4/__init__.py'
>>> 
>>> import sys
>>> sys.executable
'/server-home1/rroy/anaconda3/envs/psi4/bin/python'

Well, the good news is that when I ran this on the computer that has Jupyter (and a screen), it works fine. Downside is that that computer is running the next generation driver that’s still in branch (https://github.com/psi4/psi4/pull/1351) and a little rough for public use. So long-term this problem goes away. Short term, I still need to investigate with a current master psi. Thanks for the report.

Hi. Was this problem ever solved? I’ve just got the same error message when running the tutorial from http://www.psicode.org/psi4manual/1.2/psiapi.html on a Jupyter notebook and the 1.3.2 version (which seems the latest version from conda, psi4 channel). The error does not arises when using version 1.2.1 (or when using analytical second derivatives also with 1.3).

Thanks!

It definitely wasn’t solved by me. @loriab?

Hello.

I’ve checked that the problem seems to be that class iterator (__getitem(), for :py:class:~psi4.core.CdSalcList) are not properly handled within Jupyter (while they are correctly managed when called from python scripts). This is likely a problem on Jupyter rather than on Psi4. The workaround I’ve found (without modifying Psi4 code, so can be patched also when installed with conda) is to modify driver.py, to change the for loops where such class is involved. So for instance, instead of:

for salc in salc_list:

(where salc_list is a ~psi4.core.CdSalcList object), we can rather use:

for i in range(len(salc_list)):
    salc = salc_list[i]

Not so pythonic but works in the Jupyter notebook. The statemenst I had to update in driver.py were:

for component in salc_list[salc_index]:
for i, salc in enumerate(salc_list):
for salc in salc_list:

Hmm, it would be worthwhile for us to double-check that CdSalcList.__getitem__ in Psi’s pybind code is right before assuming Jupyter is the issue…

Thanks for the investigation! Normally I’d look into this myself, but I have a couple things to finish up before I’ll have time to look into thid.

what did you change these statements in driver.py to?

for component in salc_list[salc_index]:
for i, salc in enumerate(salc_list):
for salc in salc_list:

We have this fixed in developer Psi. We’d like to not need these hacks, but they preserve Jupyter at least.

You can see the changes we made here.

When I used the development version of psi4 downloaded using conda it still didn’t work.
Although I didn’t use the branch that loriab linked earlier in this thread

What version number does it say you have in the output file?

The version I am running is 1.4a2.dev78