Fatal Error: Wavefunction::Ca: Unable to obtain MO coefficients

1

I am trying to calculate the dipole moment and energy of a diatomic molecule, magnesium fluoride, using PsiAPI. The manual formats the code for this as follows:

E, wfn = energy('hf', return_wfn=True)
oeprop(wfn, 'DIPOLE', 'QUADRUPOLE', title='H3O+ SCF')

My code is very similar to this, it’s laid out as follows:

E, wfn = psi4.energy('scf', bsse_type='nocp', return_wfn='on')
data[1,i] = 27.211*E
data[9,i] = psi4.oeprop(wfn, 'DIPOLE')

(Note that both the energy and the dipole moments are being calculated at different molecular separations using a for loop, and the results are being saved as entries in the data matrix.) However, when I try to run this code, I get the following error message:

Fatal Error: Wavefunction::Ca: Unable to obtain MO coefficients.
Error occurred in file: D:\a\1\s\psi4\src\psi4\libmints\wavefunction.cc on line: 828

Does anyone know what this means, and what I can do about it? I’m a novice to Psi4 so I apologize if I require additional explanation to understand the issue here. I appreciate any response you can offer me.

2

Try with return_wfn=True like in the example. It needs to be a python Boolean.

3

I’ve replaced return_wfn='on' with return_wfn=True, but it continues to give the same error when I attempt to run it.

4

The following input file works perfectly fine for me.

molecule {
0 2
Mg
F 1 2 
}

set reference uhf 
_, wfn = energy('scf/cc-pvdz', return_wfn=True)
psi4.oeprop(wfn, 'DIPOLE')

This probably means that:
(a) you’re using an older version of Psi4 than I am, and you found a bug in an old version which is not present in newer versions, or
(b) there are other parts of your input file that are relevant to the problem.

Try the input file I posted. If it fails, we’re in case (a). Tell me what version of Psi4 you’re using. There are instructions here, but it’s going to be in the header of your output file. If it works, we’re in case (b). Try to simplify your input file as much as you can while still maintaining the bug, and let me know what you find.

5

The code you shared worked. According to the header of the output file, the version of Psi4 I’m using is Psi4 1.4rc2.dev1.

I was wondering if the problem had to do with the fact that the code you used was written in the Psithon style while the code I am working with is PsiAPI (being able to import Psi4 as a module is very helpful for my purposes, as I need to save my calculations in an external Excel spreadsheet intermittently and I don’t know any way to do this other than with Pandas) so I tried to convert it into PsiAPI style to see if it still ran. What I ended up using was as follows:

import psi4

psi4.geometry('''
0 2
Mg
F 1 2 
''')

psi4.set_options({'reference': 'uhf'})
_, wfn = psi4.energy('scf/cc-pvdz', return_wfn=True)
A = psi4.oeprop(wfn, 'DIPOLE')

print(A)

The code does still run, and it prints the following to the screen:

*** tstart() called on DESKTOP-78AGORJ
*** at Sat Jun 19 18:02:27 2021

   => Loading Basis Set <=

    Name: CC-PVDZ
    Role: ORBITAL
    Keyword: BASIS
    atoms 1 entry MG         line   353 file C:\Users\qweyr\psi4conda\lib\share\psi4\basis\cc-pvdz.gbs
    atoms 2 entry F          line   228 file C:\Users\qweyr\psi4conda\lib\share\psi4\basis\cc-pvdz.gbs


         ---------------------------------------------------------
                                   SCF
               by Justin Turney, Rob Parrish, Andy Simmonett
                          and Daniel G. A. Smith
                              UHF Reference
                        1 Threads,    500 MiB Core
         ---------------------------------------------------------

  ==> Geometry <==

    Molecular point group: c2v
    Full point group: C_inf_v

    Geometry (in Angstrom), charge = 0, multiplicity = 2:

       Center              X                  Y                   Z               Mass
    ------------   -----------------  -----------------  -----------------  -----------------
         MG           0.000000000000     0.000000000000    -0.883986996609    23.985041697000
         F            0.000000000000     0.000000000000     1.116013003391    18.998403162730

  Running in c2v symmetry.

  Rotational constants: A = ************  B =      0.39754  C =      0.39754 [cm^-1]
  Rotational constants: A = ************  B =  11917.93052  C =  11917.93052 [MHz]
  Nuclear repulsion =   28.575569376180002

  Charge       = 0
  Multiplicity = 2
  Electrons    = 21
  Nalpha       = 11
  Nbeta        = 10

  ==> Algorithm <==

  SCF Algorithm Type is DF.
  DIIS enabled.
  MOM disabled.
  Fractional occupation disabled.
  Guess Type is SAD.
  Energy threshold   = 1.00e-06
  Density threshold  = 1.00e-06
  Integral threshold = 1.00e-12

  ==> Primary Basis <==

  Basis Set: CC-PVDZ
    Blend: CC-PVDZ
    Number of shells: 14
    Number of basis functions: 32
    Number of Cartesian functions: 34
    Spherical Harmonics?: true
    Max angular momentum: 2

   => Loading Basis Set <=

    Name: (CC-PVDZ AUX)
    Role: JKFIT
    Keyword: DF_BASIS_SCF
    atoms 1 entry MG         line   576 file C:\Users\qweyr\psi4conda\lib\share\psi4\basis\def2-universal-jkfit.gbs
    atoms 2 entry F          line   271 file C:\Users\qweyr\psi4conda\lib\share\psi4\basis\cc-pvdz-jkfit.gbs

  ==> Integral Setup <==

  DFHelper Memory: AOs need 0.002 GiB; user supplied 0.366 GiB. Using in-core AOs.

  ==> MemDFJK: Density-Fitted J/K Matrices <==

    J tasked:                   Yes
    K tasked:                   Yes
    wK tasked:                   No
    OpenMP threads:               1
    Memory [MiB]:               375
    Algorithm:                 Core
    Schwarz Cutoff:           1E-12
    Mask sparsity (%):       0.0000
    Fitting Condition:        1E-10

   => Auxiliary Basis Set <=

  Basis Set: (CC-PVDZ AUX)
    Blend: CC-PVDZ-JKFIT + DEF2-UNIVERSAL-JKFIT
    Number of shells: 58
    Number of basis functions: 182
    Number of Cartesian functions: 218
    Spherical Harmonics?: true
    Max angular momentum: 4

  Minimum eigenvalue in the overlap matrix is 2.7340938432E-02.
  Reciprocal condition number of the overlap matrix is 1.0345261447E-02.
    Using symmetric orthogonalization.

  ==> Pre-Iterations <==

  SCF Guess: Superposition of Atomic Densities via on-the-fly atomic UHF (no occupation information).

   -------------------------
    Irrep   Nso     Nmo
   -------------------------
     A1        16      16
     A2         2       2
     B1         7       7
     B2         7       7
   -------------------------
    Total      32      32
   -------------------------

  ==> Iterations <==

                           Total Energy        Delta E     RMS |[F,P]|

   @DF-UHF iter SAD:  -298.20820771513570   -2.98208e+02   0.00000e+00 
   @DF-UHF iter   1:  -298.99448330531396   -7.86276e-01   2.39410e-02 DIIS
   @DF-UHF iter   2:  -299.04427300992290   -4.97897e-02   1.88525e-02 DIIS
   @DF-UHF iter   3:  -299.08363965317528   -3.93666e-02   1.15550e-03 DIIS
   @DF-UHF iter   4:  -299.08485346792651   -1.21381e-03   3.42887e-04 DIIS
   @DF-UHF iter   5:  -299.08498835654433   -1.34889e-04   1.16774e-04 DIIS
   @DF-UHF iter   6:  -299.08499991454511   -1.15580e-05   3.75387e-05 DIIS
   @DF-UHF iter   7:  -299.08500096590313   -1.05136e-06   1.92620e-05 DIIS
   @DF-UHF iter   8:  -299.08500127486349   -3.08960e-07   5.21122e-06 DIIS
   @DF-UHF iter   9:  -299.08500129636343   -2.14999e-08   1.55245e-06 DIIS
   @DF-UHF iter  10:  -299.08500129836165   -1.99822e-09   2.71574e-07 DIIS
  Energy and wave function converged.


  ==> Post-Iterations <==

   @Spin Contamination Metric:   7.118049029E-04
   @S^2 Expected:                7.500000000E-01
   @S^2 Observed:                7.507118049E-01
   @S   Expected:                5.000000000E-01
   @S   Observed:                5.000000000E-01

    Orbital Energies [Eh]
    ---------------------

    Alpha Occupied:

       1A1   -49.089462     2A1   -26.141550     3A1    -3.822773
       1B2    -2.339360     1B1    -2.339360     4A1    -2.337315
       5A1    -1.394906     6A1    -0.515542     2B1    -0.487894
       2B2    -0.487894     7A1    -0.300817

    Alpha Virtual:

       3B1     0.030467     3B2     0.030467     8A1     0.099460
       9A1     0.162272     4B1     0.245517     4B2     0.245517
      10A1     0.317777     1A2     0.424674    11A1     0.424674
       5B1     0.478333     5B2     0.478333    12A1     0.638346
       6B1     1.653862     6B2     1.653862    13A1     1.727912
      14A1     2.506394    15A1     4.139247     2A2     4.139247
       7B1     4.150217     7B2     4.150217    16A1     4.169999

    Beta Occupied:

       1A1   -49.084999     2A1   -26.139614     3A1    -3.812142
       1B2    -2.335514     1B1    -2.335514     4A1    -2.329933
       5A1    -1.388978     6A1    -0.496825     2B1    -0.487576
       2B2    -0.487576

    Beta Virtual:

       7A1    -0.006836     3B1     0.064717     3B2     0.064717
       8A1     0.120270     9A1     0.212482     4B1     0.281497
       4B2     0.281497    10A1     0.360608     1A2     0.451897
      11A1     0.451897     5B1     0.505784     5B2     0.505784
      12A1     0.658703     6B1     1.655202     6B2     1.655202
      13A1     1.733797    14A1     2.512249    15A1     4.139320
       2A2     4.139320     7B1     4.154637     7B2     4.154637
      16A1     4.175597

    Final Occupation by Irrep:
             A1    A2    B1    B2
    DOCC [     6,    0,    2,    2 ]
    SOCC [     1,    0,    0,    0 ]

  @DF-UHF Final Energy:  -299.08500129836165

   => Energetics <=

    Nuclear Repulsion Energy =             28.5755693761800025
    One-Electron Energy =                -477.1655848586443653
    Two-Electron Energy =                 149.5050141841027482
    Total Energy =                       -299.0850012983615898

  UHF NO Occupations:
  HONO-2 :    2 B2 1.9999529
  HONO-1 :    6 A1 1.9997442
  HONO-0 :    7 A1 1.0000000
  LUNO+0 :    8 A1 0.0002558
  LUNO+1 :    3 B2 0.0000471
  LUNO+2 :    3 B1 0.0000471
  LUNO+3 :    9 A1 0.0000024


Computation Completed


Properties will be evaluated at   0.000000,   0.000000,   0.000000 [a0]

Properties computed using the SCF density matrix

  Nuclear Dipole Moment: [e a0]
     X:     0.0000      Y:     0.0000      Z:    -1.0653

  Electronic Dipole Moment: [e a0]
     X:     0.0000      Y:     0.0000      Z:    -0.6274

  Dipole Moment: [e a0]
     X:     0.0000      Y:     0.0000      Z:    -1.6927     Total:     1.6927

  Dipole Moment: [D]
     X:     0.0000      Y:     0.0000      Z:    -4.3024     Total:     4.3024


*** tstop() called on DESKTOP-78AGORJ at Sat Jun 19 18:02:28 2021
Module time:
        user time   =        nan seconds =        nan minutes
        system time =        nan seconds =        nan minutes
        total time  =          1 seconds =       0.02 minutes
Total time:
        user time   =        nan seconds =        nan minutes
        system time =        nan seconds =        nan minutes
        total time  =          1 seconds =       0.02 minutes


Properties will be evaluated at   0.000000,   0.000000,   0.000000 [a0]
OEProp: No title given, name of density matrix used for the following properties is 'SCF alpha density'
  Nuclear Dipole Moment: [e a0]
     X:     0.0000      Y:     0.0000      Z:    -1.0653

  Electronic Dipole Moment: [e a0]
     X:     0.0000      Y:     0.0000      Z:    -0.6274

  Dipole Moment: [e a0]
     X:     0.0000      Y:     0.0000      Z:    -1.6927     Total:     1.6927

  Dipole Moment: [D]
     X:     0.0000      Y:     0.0000      Z:    -4.3024     Total:     4.3024

None

The “None” at the end seems to be what it spits out when I ask it to print the dipole moment. I’ll work on simplifying the input file to see if I can isolate the bug, but in the meantime, why is it that when I ask for it to print the dipole moment, it yields “None” instead of a number?

6

Psi always converts Psithon to PsiAPI before running any computations, so I would have been surprised if Psithon vs PsiAPI was the issue.

psi4.oeprop doesn’t return anything. For API-like access to the dipole, probably the simplest option is

wfn = psi4.properties("scf", properties=['dipole'], return_wfn=True)[1]
analytic_dipole = wfn.variable("SCF DIPOLE")

A slightly more verbose option, but closer to what you originally posted, is

wfn = energy("scf", return_wfn=True)[1]
psi4.oeprop(wfn, "DIPOLE")
quadrupole = wfn.variable("DIPOLE")

Yes, the variable is saved under slightly different names in each path. I can see this being confusing, so I’ll raise this to the other developers as something to consider standardizing.

7

I have simplified the code as much as I can while still maintaining its original functionality. I’ve also applied your suggestions so that the dipole can be stored. The code now looks like this:

import psi4
import numpy as np
import matplotlib.pyplot as plt

allVals = [1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9, 2.0, 2.1, 2.2, 2.3, 2.4, 2.5, 2.6, 2.7, 2.8, 2.9, 3.0, 3.25, 3.5, 3.75, 4.0, 4.25, 4.5, 4.75, 5.0]

nrg = np.zeros(len(allVals))
dip = np.zeros(len(allVals))

psi4.set_options({'reference': 'uhf', 'basis': 'cc-PVDZ'})

for R in allVals:
    molecule = '''
        Mg 
    --
        F 1 ''' + str(R) + '''
    '''

    psi4.geometry(molecule)
    E, wfn = psi4.energy('scf', bsse_type='cp', return_wfn=True)
    nrg[allVals.index(R)] = E
    psi4.oeprop(wfn, 'DIPOLE')
    dip[allVals.index(R)] = wfn.variable('DIPOLE')
plt.plot(allVals, nrg, allVals, dip)
plt.show()

It still returns the error I mentioned: Fatal Error: Wavefunction::Ca: Unable to obtain MO coefficients. Error occurred in file: D:\a\1\s\psi4\src\psi4\libmints\wavefunction.cc on line: 828. This error occurs at line 22, psi4.oeprop(wfn, 'DIPOLE').

Out of curiosity, I tried modifying the code so that the input to the psi4.geometry() function looked more like the specifications of the molecule in your example:

molecule = '''
        0 2
        Mg 
        F 1 ''' + str(R) + '''
    '''

This gives a different error: N-Body requires active molecule to have more than 1 fragment.

On the advice of a colleague, I tried re-introducing the – between the two atoms but leaving the numbers at the beginning, as follows:

    molecule = '''
        0 2
        Mg 
        --
        F 1 ''' + str(R) + '''
    '''

This gives yet another error: Inconsistent or unspecified chg/mult: sys chg: None, frag chg: [0.0, None], sys mult: None, frag mult: [2, None].

I do not know what either of these errors mean, but the first occurs when the code tries to execute E, wfn = psi4.energy('scf', bsse_type='cp', return_wfn=True) and the second occurs one step earlier when the code tries to execute psi4.geometry(molecule). Since both of these occur before the problem line I do not know if the error would still occur under these conditions.

8

The problem is the BSSE code isn’t setting the orbitals on the wavefunction. This may seem strange, but there are all sorts of annoying definitional questions that appear if you ask how you define BSSE vs non-BSSE orbitals and dipoles. I’ll get in touch with the developers who work with BSSE more.

The “N-Body requires active molecule” error means that Psi doesn’t know how to do a BSSE computation for a system of one fragment, and the hyphen you removed was repsonsible for telling Psi there were two fragments.

In your last example, the 0 2 needs to be on the fluorine fragment, because it’s the fluorine that’s the doublet.

1 Like
9

Thank you, I appreciate that. If I can figure out how to address the BSSE problem I’ll post an update about it.

10

The issue is that with the N-body calculator is that you don’t get the actual wavefunction object(s) back, but an empty dummy.

A solution is to do the interaction energy calculation manually using this function:
https://psicode.org/psi4manual/master/api/psi4.core.Molecule.html#psi4.core.Molecule.extract_subsets

If your molecule object is called mol then this should work if I didn’t introduce typos.

mAcp =mol.extract_subsets(1,2)
eA = psi4.energy(...,molecule=mAcp)
mBcp =mol.extract_subsets(2,1)
eB = psi4.energy(...,molecule=mBcp)
edimer, wfn  = psi4.energy(...,return_wfn=True, molecule=mol)
eCP = edimer - eA - eB
analytic_dipole = wfn.variable("SCF DIPOLE") # dipole is calculated by default, no need to call oeprop
print(f" E(int,cp) = {eCP}")
11

An important warning if you go that route: the dipole you’re computing is not BSSE-corrected. If you want a BSSE corrected dipole, I think it’s doable, but the code is going to be clunky.

12

This comment is for the general case of N-body property calculations (many body expansion).

For your dimer, you don’t need to worry about BSSE, it is the full cluster in the MBE sense.

13

Finally got a chance to try running this. It gives an AttributeError. AttributeError: 'psi4.core.Molecule' object has no attribute 'extract_subset' I’ve defined the script I’m using to test this very similarly to yours, for clarity I’ll include the whole thing as it isn’t very long:

import psi4
import numpy as np
import matplotlib.pyplot as plt

allVals = [1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9, 2.0, 2.1, 2.2, 2.3, 2.4, 2.5, 2.6, 2.7, 2.8, 2.9, 3.0, 3.25, 3.5, 3.75, 4.0, 4.25, 4.5, 4.75, 5.0]

nrg = np.zeros(len(allVals))
dip = np.zeros(len(allVals))

psi4.set_options({'reference': 'uhf', 'basis': 'cc-PVDZ'})

for R in allVals:
    molecule = '''
        Mg 
        --
        0 2
        F 1 ''' + str(R) + '''
    '''

    mol = psi4.geometry(molecule)
    mAcp = mol.extract_subset(1,2)
    eA = psi4.energy('scf', molecule=mAcp)
    mBcp = mol.extract_subset(2,1)
    eB = psi4.energy('scf', molecule=mBcp)
    edimer, wfn = psi4.energy('scf', return_wfn=True, molecule=mol)
    nrg[allVals.index(R)] = edimer - eA - eB
    dip[allVals.index(R)] = wfn.variable("DIPOLE")

plt.plot(allVals, nrg, allVals, dip)
plt.show()
14

The method name is extract_subsets, with the s on the end. (It’s plural because you can pass lists to get out multiple fragments.)

15

I’ve added an s to the end of both extract_subset calls (so they now say extract_subsets). I’m now getting a different error:

Exception has occurred: KeyError
'psi4.core.Wavefunction.variable: Requested variable DIPOLE was not set!\n'

Do you know what this means? It just seems like I’m running into one issue after another here.

16

This particular error has to do with something I said previously: whether the variable is saved as “DIPOLE” or “SCF DIPOLE” depends on how you pass it. I’ll make an issue about standardization, because it shouldn’t be this confusing.

After changing the variable to SCF DIPOLE, there’s still an error: SCF DIPOLE has three components, so you either need to pick out the one non-zero z-component or change the dimensions of dipole to len(allVals) x 3. Based on your matplotlib, I’m assuming you want the z-component. To fix both of the errors so far, change wfn.variable("DIPOLE") to wfn.variable("SCF DIPOLE")[2].

That said, there’s still another error: SCF convergence difficulties starting at 2.8 A. While I’m not surprised to see convergence issues along bond dissociation, I am surprised that my normal tricks aren’t solving them. I can investigate this, but don’t expect results for a few more days.

17

Thank you very much, with that suggestion the simplified version of my code now runs to completion! I suspect that the full version will still encounter some rather pesky errors but I’ve already opened up a thread for one of those.

I think I’ve actually run into those convergence problems already, although it happened when I was trying to calculate the energy, not the dipole moment. I’ve been trying to temporarily deal with that using an except statement to just define an arbitrary value for the energy so it can move on to the next step whenever it fails to converge instead of causing the whole script to fail; I was hoping to try to figure out that issue later once some of these more fundamental problems are fixed. Still, it’s refreshing to hear that I’m not the only one encountering this problem. If you are able to devise a workaround, please let me know. I don’t mind waiting.

Thank you all for your help, I really appreciate it! You’ve allowed me to make a big step in ironing out the problems in the code today.

18

Out of curiosity, is there such a thing as a CCSD(T) dipole moment calculation, or does the CCSD(T) method only improve on SCF when we’re calculating energy? Elsewhere in my code I’m calculating energies using the CCSD(T) method and I’m wondering if it is meaningful to attempt to pull out the dipole moment here.

19

Available properties are listed here:
https://psicode.org/psi4manual/master/prop.html

1 Like
20

Optimizing orbitals at some of these more stretched geometries is a difficult numerical problem (the technobabble is that the molecular orbital hessian is near-singular and has a negative eigenvalue near the minimum). Psi doesn’t have the technology to optimize orbitals reliably there.