15
$\begingroup$

I'm doing single point energy calculations of FeS using ORCA. I originally used DFT with a variety of functionals, and these calculations all took less than a minute. Now I'm attempting to run CASSCF calculations using 8 active orbitals, 12 active electrons, and the aug-cc-pVTZ basis set. It's been running for 3 days. Is this typical, or should I be attempting to change tolerances and other settings? How long would one anticipate something like this to run? In addition, how can I (generally speaking) evaluate how close a given calculation is to finishing?

$\endgroup$
  • 1
    $\begingroup$ Would you consider remove diffuse function, which may save some time. In addition, molpro's casscf is typically faster than orca. Though dmrg can be faster, I heard it only becomes useful for beyond ~(12,12) $\endgroup$ – Rodriguez Jun 23 '16 at 16:13
  • $\begingroup$ I can definitely try that! I'm replicating data from a paper, so it'll just be a matter of seeing how close it gets. Removing the diffuse function would be making my basis set cc-pVTZ instead of aug-cc-pVTZ, correct? And that's good to know about Molpro, although since I just started and am starting with ORCA, I'm going to try to stick to that for the time being. $\endgroup$ – Argon Jun 23 '16 at 16:22
  • $\begingroup$ Yes. I mean cc-pVTZ.. $\endgroup$ – Rodriguez Jun 23 '16 at 16:25
  • 1
    $\begingroup$ Did you check the partial values of the ocupation numbers? Close to full or empty ocupation can give convergence problems $\endgroup$ – user1420303 Jun 23 '16 at 23:41
  • $\begingroup$ No, I haven't. I'll have to look into how to do that because I'm not sure how yet, but I'll give that a go. Thank you! $\endgroup$ – Argon Jun 23 '16 at 23:46
24
$\begingroup$

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.

orca_casscf.txt energies

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:

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

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

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

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

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

$\endgroup$
  • $\begingroup$ I can't thank you enough for such a thorough response. I am taking my time to digest everything you've said, and it's helping a lot, so I really appreciate it. I realize in the paper that they used a (12,14) active space. I thought it would be less costly to use a smaller active space, but that clearly wasn't the case. For some reason (12,10) ran extremely quickly and (12,14) is still running after about a week. Still, perhaps that means the orbitals chosen as the active orbitals by ORCA aren't correct. I will be sure to investigate further. Thanks for getting me started! $\endgroup$ – Argon Jul 7 '16 at 4:55
  • 1
    $\begingroup$ What I tried to convey is that everything is a balance. The (12,8) didn't run because there was no way to optimize the wavefunction. The (12,10) ran in reasonable time, but you could probably get a better energy from rotating your orbitals. The (12,14) is going to take a while (even with good orbitals) because the number of CSFs you're generating is huge, and this will require a lot of memory, plus your comment about the identity of the orbitals. $\endgroup$ – pentavalentcarbon Jul 7 '16 at 15:57
  • $\begingroup$ @AndrewR. Also, feel free to ping me in chat if you have any questions. $\endgroup$ – pentavalentcarbon Jul 9 '16 at 23:34
  • $\begingroup$ thank you for the help again. I've taken your comment into consideration and have a few (hopefully short) follow-up questions. I'm new to Stack Exchange so don't know how to ping you in chat. I created a room, but I'm not sure I was supposed to do that... $\endgroup$ – Argon Jul 12 '16 at 0:01
  • $\begingroup$ It's fine; to ping me in chat, go here (chat.stackexchange.com/rooms/3229/the-periodic-table), type '@pentavalentcarbon' and then some text, and I'll get a notification. I can't find the room you created on the chat page; if you'd rather chat there, you'll have to send me a link. $\endgroup$ – pentavalentcarbon Jul 12 '16 at 3:50
3
$\begingroup$

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.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.