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

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.

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

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

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.

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

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

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 <==

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?

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")


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.

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.

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

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

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}")


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.

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.

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()


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.)

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.

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.

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.

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.

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

1 Like

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.