Restarting an FCI calculation in PySCF

Question

My question is, if there is a way to continue an earlier FCI calculation in PySCF.

Imagine you initially just need energies converged up to a certain threshold:

  cisolver = pyscf.fci.FCI(an_scf_object)
  cisolver.conv_tol = 1e-5
  e, fcivec = cisolver.kernel(verbose=5)

and later you decide that you need very accurate energies and CI vectors:

  cisolver = pyscf.fci.FCI(an_scf_object)
  cisolver.conv_tol = 1e-9
  e, fcivec = cisolver.kernel(verbose=5)

What properties could you write down after the initial run in order to restart from there?

Answer 1

What you can do, is save the CI vector of your calculation and then use it as starting point in a subsequent calculation, like so:

e, c = cisolver.kernel(verbose=5,ci0=c0) 

where c0 is a numpy array containing the solution of another run. I tested this for carbon in a STO-3G basis. First ran to a convergence of 1e-5, saved the vector, then converged with saved vector to 1e-9. This I compared to a run from scratch which was converged to 1e-9. The continuation took 4.15s to converge while the calculation from scratch took 5.21s.

Answer 2

Preliminaries

The hint by Sverus Snape to use a numpy array was helpful, but how does one do that?

I found the answer at How to save and load numpy.array() data properly? and since there's also so many different options given there (e.g. np.fromfile, np.loadtxt, np.load, etc.) I would like to provide people with the solution that worked for me, and a minimum working example.

Benchmark calculation (without restarting)

A calculation of the CISD energies for the lowest three electronic states of Ne in cc-pV5Z using this input with conv_tol=1e-9 required 9 Davidson iterations to converge the ground state energy to this precision, and more than 30 Davidson iterations to converge the third highest of these electronic states:

davidson 34 5  |r|= 6.05e-05  e= [-0.33719767  0.80457309  0.80457311]  max|de|= -2.83e-09  lindep= 0.913
davidson 35 6  |r|= 3.02e-05  e= [-0.33719767  0.80457309  0.80457311]  max|de|= -1e-09  lindep= 0.738
root 2 converged  |r|= 2.1e-05  e= 0.8045731078550442  max|de|= -1.36e-10
converged 36 7  |r|= 2.1e-05  e= [-0.33719767  0.80457309  0.80457311]  max|de|= -1.36e-10
RCISD converged
RCISD root 0  E = -128.8839677954546
RCISD root 1  E = -127.7421970392584
RCISD root 2  E = -127.7421970216469

Partial calculation (saves CI vector for later)

The example below does the same calculation but wtih conv_tol=1e-5 and it saves the numpy array of the CI vector. I ended up using the answer with the second highest net score on the above StackOverflow thread, because it was more straightforward, and because binary .npy files are more practically sized than ASCII .txt files (a solution with .h5 would be even better, and for the PySCF coupled cluster code, they do document such a solution). The command that I used for saving the CI vector is in the last line below:

#!/usr/bin/env python

import numpy as np
import pyscf
from pyscf import gto, scf, ao2mo, fci,ci

mol = pyscf.M(atom = 'Ne 0 0 0',basis = 'cc-pv5z',verbose=5,output='out_1e-5_direct.txt')
mhf = scf.RHF(mol).run()
mci = ci.CISD(mhf).set(conv_tol=1e-5,nroots=3)
e, civec = mci.kernel()
np.save('civec.npy', civec)

Because of the lower convergence criterion, the ground state energy converged after only 5 Davidson iterations, and the highest of the three states sought had its energy converge at around the 19th Davidson iteration (all three energies are also slightly higher than in the previous calculation):

davidson 17 6  |r|= 0.00551  e= [-0.33719763  0.80457337  0.80459828]  max|de|= -1.57e-05  lindep= 0.691
davidson 18 7  |r|= 0.00339  e= [-0.33719763  0.80457337  0.80459439]  max|de|= -3.89e-06  lindep= 0.443
root 2 converged  |r|= 0.00128  e= 0.8045935429617836  max|de|= -8.5e-07
converged 19 8  |r|= 0.00132  e= [-0.33719763  0.80457337  0.80459354]  max|de|= -8.5e-07
RCISD converged
RCISD root 0  E = -128.8839677596848
RCISD root 1  E = -127.7421967633118
RCISD root 2  E = -127.7421765865402

My directory then had the following checkpoint file (not human-readible, since it's a binary file):

-rw-r----- 1 nike nike 4448072 Jan  4 17:42 civec.npy

Calculation that reuses the saved CI vector

Now I'll recover the "benchmark" energies (conv_tol=1e-9) using the .npy file that was saved in the last line of the script in the previous section; it is reused in the last line of the following script:

#!/usr/bin/env python
# Author: Nike Dattani, [email protected]

import numpy as np
import pyscf
from pyscf import gto, scf, ao2mo, fci,ci

mol = pyscf.M(atom = 'Ne 0 0 0',basis = 'cc-pv5z',verbose=5,output='out_1e-9_restarted.txt')
mhf = scf.RHF(mol).run()
mci = ci.CISD(mhf).set(conv_tol=1e-9,nroots=3)
e, civec = mci.kernel(ci0=np.load('civec.npy'))

This gives roughly the same energies from the "benchmark" calculation, but with only about 2 Davidson iterations for the ground state energy, and about 12 iterations for the third energy sought:

davidson 0 3  |r|= 0.00132  e= [-0.33719763  0.80457337  0.80459354]  max|de|= 0.805  lindep= 0.816
davidson 1 6  |r|= 0.000962  e= [-0.33719767  0.80457312  0.80459334]  max|de|= -2.48e-07  lindep= 0.506
root 0 converged  |r|= 2.04e-05  e= -0.33719766591648676  max|de|= -8.9e-10
davidson 2 9  |r|= 0.000859  e= [-0.33719767  0.80457309  0.80459321]  max|de|= -1.32e-07  lindep= 0.63
davidson 3 11  |r|= 0.00165  e= [-0.33719767  0.80457309  0.80459283]  max|de|= -3.75e-07  lindep= 0.796
root 1 converged  |r|= 5.81e-06  e= 0.8045730910409581  max|de|= -6.13e-11
davidson 4 13  |r|= 0.00418  e= [-0.33719767  0.80457309  0.80458896]  max|de|= -3.88e-06  lindep= 0.939
davidson 5 14  |r|= 0.00519  e= [-0.33719767  0.80457309  0.80457749]  max|de|= -1.15e-05  lindep= 0.892
davidson 6 15  |r|= 0.00117  e= [-0.33719767  0.80457309  0.80457335]  max|de|= -4.14e-06  lindep= 0.914
davidson 7 16  |r|= 0.00041  e= [-0.33719767  0.80457309  0.80457318]  max|de|= -1.66e-07  lindep= 0.775
davidson 8 3  |r|= 0.00041  e= [-0.33719767  0.80457309  0.80457318]  max|de|= 9.11e-11  lindep= 0.999
davidson 9 4  |r|= 0.000371  e= [-0.33719767  0.80457309  0.80457314]  max|de|= -3.54e-08  lindep= 0.692
davidson 10 5  |r|= 0.000219  e= [-0.33719767  0.80457309  0.8045731 ]  max|de|= -4.39e-08  lindep= 0.929
davidson 11 6  |r|= 7.38e-05  e= [-0.33719767  0.80457309  0.80457309]  max|de|= -7.99e-09  lindep= 0.781
root 2 converged  |r|= 3.04e-05  e= 0.8045730916435243  max|de|= -7.2e-10
converged 12 7  |r|= 3.04e-05  e= [-0.33719767  0.80457309  0.80457309]  max|de|= -7.2e-10
RCISD converged
RCISD root 0  E = -128.8839677954263
RCISD root 1  E = -127.7421970385624
RCISD root 2  E = -127.7421970378585

Extra remarks

It's remarkable that when following the "save and restart" method, the total number of Davidson iterations is actually smaller than what was required when setting conv_tol=1e-9 directly in one calculation.

This procedure might also be useful for finding higher excited state energies if your nroots parameter was set undesirably low during your first calculation.

Availability of above input and output files

The input and output files for all three calculations in this answer are provided here. The files are described below:

Calculation	Input file	Output file
Benchmark calculation (conv_tol=1-e9) with no restart	inp_1e-9_direct.py	out_1e-9_direct.py
Partial calculation (conv_tol=1-e5), saves the checkpoint file	inp_1e-5_direct.py	out_1e-5_direct.py
Benchmark recovery (conv_tol=1-e9), using checkpoint file	inp_1e-9_restarted.py	out_1e-9_restarted.py