Retry un-converged SCF from Python

I am interested in the use case described in this old post: Save orbitals and restart SCF
Namely, restarting a failed SCF job, with the old calculation as the wavefunction guess, but adding slow-and-steady convergence options such as SOSCF or damping.

The code included in that post (translated to Python) does not work for me; in fact, GUESS: READ seems to do nothing when a restart is initiated in this way. The SCF restart happens, but I just get SCF Guess: Superposition of Atomic Densities via on-the-fly atomic UHF (no occupation information) as usual, without even an error message saying there is no wavefunction file (there is, I checked using os.path.exists()).

This example code should be complete and runnable. Note that, when this code is run, the second SCF run in the resulting file does not use the guess from the first. SOSCF does kick in, as desired, but it seems quite wasteful not to use the existing wavefunction guess.

'''
Test psi4 SCF restart capabilities
'''
import os
import time

import numpy as np
import psi4

def run_psi4(name, atomic_numbers, xyz, net_charge=0, theory='wB97X-D3BJ/def2-TZVPD',
             ram_gb=4.0):
    '''
    Run psi4
    '''
    outfile = name + '.psi4'
    psi4.core.set_output_file(outfile, False)
    atomic_numbers = np.array(atomic_numbers)
    n_electrons = np.sum(atomic_numbers) - net_charge
    psi4.set_memory('%f GiB' % ram_gb)
    geom_string = '\n'.join(['%d %.8f %.8f %.8f' % (e, x, y, z) for e, (x, y, z) in zip(atomic_numbers, xyz)])
    geom_string = str(net_charge) + ' 1\n' + geom_string
    mol = psi4.geometry(geom_string)  # in angstroms by default
    settings = {
        'scf_type': 'DF',  # density fitting, may use RAM or disk (needs fast disk)
        'maxiter': 10,  # for testing, make sure we don't converge first time
    }
    psi4.set_options(settings)
    start_time = time.time()
    try:  # allow SCF failure without crashing
        grad, wfn = psi4.gradient(theory, return_wfn=True)  # in Hartree/Bohr
    except psi4.driver.SCFConvergenceError as ex:
        print(repr(ex), '- retrying', flush=True)
        # try all obvious ways of writing the wfn file
        wfn_file = ex.wfn.get_scratch_filename(180)
        ex.wfn.to_file(wfn_file)
        wfn_file = ex.wfn.get_scratch_filename(180) + '.npy'
        ex.wfn.to_file(wfn_file)
        print(f'{wfn_file = }, {os.path.exists(wfn_file) = }', flush=True)
        psi4.set_options({
            'SOSCF': True,  # 2-3x slower, but steadier convergence
            'guess': 'read',
            'maxiter': 100,
        })
        try:  # try again with new options
            grad = psi4.gradient(theory)  # in Hartree/Bohr
        except psi4.driver.SCFConvergenceError as ex:
            psi4.core.clean()
            psi4.core.clean_options()
            print(repr(ex), '- stopping', flush=True)
            return  # skip this system
    psi4.core.print_variables()  # list variables explicitly in output file
    psi4.core.clean()  # clean up old data & files for next round
    psi4.core.clean_options()
    elapsed_time = time.time() - start_time
    print(name, n_electrons, 'electrons', elapsed_time, 'seconds', flush=True)

if __name__ == '__main__':
    n_cores = 8
    ram_gb = 32.0
    psi4.set_num_threads(n_cores)
    name = 'Li-metal-cluster'
    atomic_numbers = np.array([3] * 16)
    xyz = np.array([[ 4.00,    0.00,    0.25],
                    [-2.00,   -0.00,   -1.75],
                    [-2.00,   -0.00,    2.25],
                    [-4.00,   -0.00,    0.25],
                    [-2.00,    2.00,    0.25],
                    [-2.00,   -2.00,    0.25],
                    [ 2.00,    0.00,   -1.75],
                    [ 0.00,    0.00,   -3.75],
                    [-0.00,    2.00,   -1.75],
                    [ 0.00,   -2.00,   -1.75],
                    [ 2.00,    0.00,    2.25],
                    [ 0.00,    0.00,    0.25],
                    [-0.00,    2.00,    2.25],
                    [ 2.00,    2.00,    0.25],
                    [ 0.00,   -2.00,    2.25],
                    [ 2.00,   -2.00,    0.25]])
    net_charge = 0
    run_psi4(name, atomic_numbers, xyz, net_charge,
                 'wB97X-D3/cc-pVTZ', ram_gb)

To be honest, I have been struggling with the Psi4 restart functionality since I first started using Psi4 a couple of years ago. As far as I have searched, documentation for the 180 file and wavefunction guesses is still not complete, despite similar forum posts like Read in from scf with psi4.core.Wavefunction.from_file()? - #2 by jmisiewicz.

  • No, the orbital guess is still not properly documented. While I’m unhappy about it, the frustrating reality is that the developers’ time is so limited that the sorely needed documentation never becomes important enough to do now. I personally am facing the choice between working on documentation and having the flexibility to respond to bug reports (like this) that don’t clearly belong to any other developer. Clearly, the correct choice is to bother me about orbital reads until it’s better for me to write the documentation than it is to keep answering question.

  • You’ve stumbled across a new bug, thanks for that. I’ve filed it here: `guess` Not Updated After Failed Job · Issue #2458 · psi4/psi4 · GitHub. This is an “I’ll deal with it over the weekend, or earlier if my schedule unexpectedly lightens” severity bug.

  • 1.6 (targeted for release in May) adds ADIIS and EDIIS support to Psi, which should hopefully prevent the SCF from failing in the first place.

For energy calculation the orbital guess has been simplified as described here:HF: Hartree–Fock Theory

However this does not work for gradient calculations.

I’m glad to hear you’re on top of it. I wasn’t sure whether this guess behavior was a bug or intended-but-undocumented. Thanks for the github issue, and for bringing up the new SCF convergence features expected in v1.6. I’m running Psi4 on large numbers of metal clusters and seeing about 35% fail in SCF convergence, as compared to negligible failure rates with organic systems, so convergence is on my mind right now.

below is a quick example how to use the write_orbitals functions and how to rename and copy the numpy file so that psi4 sees it.

To get change the SCF guess setting to change after the not-converged-exeception it appears one needs to change the local option not the global one.

orbital_file='mos' # written where psi4 is executed.
import shutil

try:
  E, wfn = gradient('hf/sto-3g',write_orbitals=orbital_file, return_wfn=True)
except psi4.driver.SCFConvergenceError as ex:
  psi4.core.set_local_option('SCF', 'GUESS' ,'READ')
  scr = psi4.core.IOManager.shared_object().get_default_path() # psi4 scratch dir
  molecule = psi4.core.get_active_molecule()
  molname =  os.path.split(os.path.abspath(psi4.core.get_writer_file_prefix(molecule.name())))[1]
  target = os.path.join(scr, molname + ".180.npy") # rename
  origin = orbital_file + ".npy"
  shutil.copy(origin, target)  # copy
  psi4.core.print_out(f" \n Copying restart file <{origin}> to <{target}> \n")
  E, wfn = gradient('hf/sto-3g', return_wfn=True)

In your script, maybe using psi4.core.set_local_option('SCF', 'GUESS' ,'READ') is enough.

@hokru yes, I can confirm that adding
psi4.core.set_local_option('SCF', 'GUESS', 'READ')
after my psi4.set_options call works around this bug and makes the guess load correctly in my code example above.

Very interesting about using restart_file, thank you. I saw that documentation and tried using restart_file with gradient(), not knowing that it only works with energy(). Seeing it not work, I looked around and got the impression from a forum post here (Restarting a SCF with another reference SCF wavefunction - #7 by jmisiewicz) that restart_file was an older mechanism, rather than preferred.

At that time I wrote that message, it was. :slight_smile: restart_file was only made operational again in July 2021, thanks to @hokru. I completely forgot about that.