Chemistry - What are typical runtimes for CASSCF calculations?
Solution 1:
There is nothing trivial about MCSCF calculations because it is hard to predict a priori how long a calculation will take. There are well-defined equations for calculating how many determinants
$$ D(n,N,S) = \binom{n}{N/2+S} \binom{n}{N/2-S} $$
or configuration state functions (CSFs)
$$ D(n,N,S) = \frac{2S+1}{n+1} \binom{n+1}{N/2-S} \binom{n+1}{N/2+S+1} $$
will be required for your problem; $n$ is the number of active orbitals, $N$ is the number of active electrons, and $S$ the total spin. The calculation cost will scale with this number and can be further reduced with the use of point group symmetry. (Actually, depending on your system and the property of interest, calculations without symmetry can be meaningless, like for excited states and their energies.)
These equations should make it clear that the cost grows very rapidly, considering that for CASSCF, a full CI is being performed within your chosen active space. This is not what makes MCSCF calculations tricky in practice. Here is a plot of the total energy after the first 30 iterations.
This is because the occupation numbers are very close to 0/1/2, and they aren't changing with each iteration. Notice that the occupation numbers have inverted themselves, and repeated at every step is Rot=23,19. Orbital 19 is within the active space, but 23 is external. This is an indication that the orbitals in the active space are not the correct ones. From the paper:
the 12 or 13 valence electrons of the neutral or anion, respectively, were distributed among 14 valence active orbitals: 3p of sulphur plus 3d, 3d' and 4s orbitals of iron. With this active space we found that all important correlating orbitals are included in the active space.
MCSCF is not a black box; it's one of the few quantum chemical methods left that requires chemical intuition to even get started. This a (12,14) active space; the one in your calculations is (12,8). The procedure I use for running MCSCF calculations is generally:
Converge a single-reference calculation. This could be HF, DFT, MP2, etc. The usual recommendation is to form MP2 natural orbitals, then check the occupations numbers. If the occupation deviates more than 0.02 from an integer value, then it should probably be in the active space. This is in addition to considering the electrons and orbitals you wanted to correlate to begin with, along with antibonding partners. For example, in ethylene this is probably the pi bonding/antibonding pair, leading to a (2,2), and for 3d metal complexes this is whatever coordination (anti)bonds are present, with potentially the double d shell. I've had some luck with using DFT orbitals for metal complexes, where perturbation theory might be a bad idea due to near-degeneracies or cost. The goals here are to confirm your active space and provide the best possible starting orbitals for MCSCF.
Look at your orbitals. Of course, read the Mulliken or Lowdin population analysis, but from personal experience, I've caught many potential mistakes by doing this. It also helps when visualizing how the single-state, single-reference orbitals transform to the natural orbitals (potentially state-average) that result from MCSCF calculations.
Reorder your orbitals if necessary before running the MCSCF. Except for the most trivial cases, the orbitals that belong in the active space are not the orbitals that come straight out of your single-reference calculation. If they were, then you might not be running MCSCF to begin with. I also tend to look at the reordered orbitals one more time before doing MCSCF to confirm I haven't made any mistakes.
Read in the previous set of (reordered) orbitals into a separate MCSCF calculation and cross your fingers. Monitoring the occupation numbers of your active space is a good guide for seeing how the calculation is converging.
Do whatever post-processing you need for your project. For me, this is looking at final occupation numbers, the orbitals themselves, and then doing MR-CISD or some other multireference correlated calculation.
Doing an auto-occupation procedure without starting from previously-converged orbitals means that 1. you don't have any guarantees about the contents of your active space, 2. you start from orbitals that are close to Hartree-Fock in quality, and 3. the calculation has to do far more work in optimizing both sets of variational parameters (the MO coefficients $\{C\}$ and the CI coefficients $\{c\}$) since the internal and external orbitals (where the CI coefficients are frozen) are still far from convergence using a single-reference method.
To walk through part of this workflow, I ran a PBE0/aug-cc-pVTZ calculation (converges much faster than the SCF before MP2) to look at the Lowdin MO populations:
! pbe0 aug-cc-pvtz cc-pvtz/jk ri rijk tightscf usesym zora
%output
print[p_orbpopmo_l] 1
end
* xyz 0 5
Fe 0.000000 0.000000 0.000000
S 1.960000 0.000000 0.000000
*
The MO indices are S 3s: 14, Fe 4s: 20, Fe 3d: 15,16,17,18,19, and Fe 3d': 33,34,36-40. Automatic occupation would have gotten the occupieds right, but the virtuals wrong. I'm not sure which of the 7 virtual MOs that might be the Fe 3d' are the correct ones; there is probably some trial-and-error here, but these should be visualized. Another point is that I've only looked at the ordering for alpha-spin orbitals; for difficult cases such as transition metals, the spatial ordering for beta-spin orbitals can be very different. The reordering for alpha- and beta-spin orbitals might be different. Again, this is speaking from experience.
The final exhibit is the RI-MP2/aug-cc-pVTZ natural orbitals that might enter an MCSCF calculation based on their occupation numbers:
N[ 15]( A2) = 1.97761944
N[ 16]( A1) = 1.97753198
N[ 17]( B2) = 1.97009578
N[ 18]( B1) = 1.97009578
N[ 19]( A1) = 1.03828526
N[ 20]( A1) = 0.99215102
N[ 21]( B1) = 0.98317529
N[ 22]( B2) = 0.98317529
N[ 23]( A2) = 0.02013554
N[ 24]( B2) = 0.01947713
N[ 25]( B1) = 0.01947713
N[ 26]( A1) = 0.01816184
To be thorough, this should be combined with looking at the Lowdin MO populations of these, which might require reading them into another (SCF) calculation and not performing any iterations.
However, it is well-known that converging MCSCF equations is generally difficult. In my limited experience, DALTON, GAMESS, and Molcas have less trouble than ORCA. I haven't used Molpro for MCSCF. Which program you choose dictates what post-MCSCF calculations can be performed. The DIIS algorithm in ORCA, in particular, can stall even when it shouldn't. If you can afford it, the Newton-Raphson algorithm (switchstep nr), where the active space orbital Hessian is calculated at each step, works very well to force convergence once SuperCI is done. ORCA makes it easy to get started with CASSCF, and it has plenty of knobs, but it has a very slow integral engine, doesn't prevent you from doing stupid things like GAMESS, and more advanced features are poorly documented. The workflow described above (namely the plotting orbitals and rotating them) becomes very tedious as well, but that has more to do with MCSCF in general.
Here are some good resources that describe both practical and theoretical aspects of performing MCSCF calculations: 1, 2, 3, 4.
Solution 2:
It seems increasing the active space helped make the CASSCF calculation converge. I've had success with increasing the number of active orbitals from 8 to 10. It only took an hour or so.