Setting thresholds for integral in CASSCF/RASSCF

[ankit7540]

Dear all,
I am running MCSCF calculations for finally obtaining good electronic wavefunctions for later response calculations. Convergence of the electronic wavefuction is slow and I am looking for ways to speed things up. I wonder if reducing the threshold for the integrals maybe a good way.

From the manual of DALTON,

Code: Select all

*SUPINT


.THRESH
READ (LUCMD,*) THRSUP
Threshold for the supermatrix integrals. Read one more line containing the threshold.
Default is the same as the threshold for discarding the two-electron integrals (see the
chapter describing the molecule input format, Ch. 26).

is one way to set the threshold. Next, in Ch 26 says that in the molecule input file we can define the threshold as

Code: Select all

Integrals (Real). Indicates the threshold for which integrals smaller than this will
be considered to be zero. If not specified, a threshold of 1.0D-15 will be used. A
threshold of 1.0D-15 will give integrals correct to approximately 1.0D-13.

My questions.

• 1. Are both parameters need to defined during the calculation or one is enough.

• 2. How to properly give input in the *.mol file.

Why I am eager to try this.
For O2 I received this output (full file attached), at the end of the electronic part of the calculation,

Code: Select all

 Magnitude of CI coefficients
 ============================
 ( Ranges are relative to norm of vector :  1.00E+00 )

  10- 1 to 10- 0         6    0.88051459E+00      0.8805145941
  10- 2 to 10- 1       152    0.87662190E-01      0.9681767840
  10- 3 to 10- 2      3156    0.27643054E-01      0.9958198378
  10- 4 to 10- 3     42066    0.37854091E-02      0.9996052469
  10- 5 to 10- 4    292457    0.38256489E-03      0.9999878118
  10- 6 to 10- 5    609481    0.12082055E-04      0.9999998939
  10- 7 to 10- 6    357798    0.10591637E-06      0.9999999998
  10- 8 to 10- 7     64140    0.23031285E-09      1.0000000000
  10- 9 to 10- 8      8060    0.26174301E-12      1.0000000000
  10-10 to 10- 9      4774    0.95959001E-15      1.0000000000
  10-11 to 10-10      4490    0.10694157E-16      1.0000000000
  Number of coefficients less than 10^-11 times norm is                   3198

  Total CPU  time used in SIRIUS :  15 hours 55 minutes 26 seconds
  Total wall time used in SIRIUS :   9 hours 24 minutes 15 seconds

Very small coefficients (<1e-9) maybe be neglected. That's my guess but I am seeking suggestions from the developers. Thank you.

[ankit7540]

Update 1 : After some testing,

Code: Select all

*SUPINT
.THRESH

does not significantly affect run times. Also it has no effect on the CI coefs.

Update 2: After some testing,

Code: Select all

.Integrals

in the *.mol files also does not significantly reduce run times.

[taylor]

It would of value to see exactly what calculation you are trying to run on O2 that converges slowly and that takes a very substantial amount of time (for a diatomic). It might be possible to improve on this.

I do not see that reducing integral accuracy (which is not a very good idea for reasons mentioned below) can improve convergence, and I do not see that it will reduce the computer time per iteration. In a given iteration the integrals over symmetry orbitals are transformed to the current MO basis, and then used to solve the MCSCF equations to obtain improved MOs and CI coefficients. Neglecting small integrals does not reduce the time required for modern symmetry-blocked integral transformations, and of course once the MO integrals are available the computer time required will be completely independent of any neglect of integrals in the symmetry orbital basis.

Why did I say reducing integral accuracy was not a very good idea above? If the accuracy is not reduced by much, it probably has little effect for the calculation of MCSCF energies and first-order properties. It may have a larger effect on higher-order/response properties. If the accuracy is significantly reduced, however, all sorts of problems may appear, including variational collapse because the electron repulsion contribution is not given correctly.

The number of configurations in your calculation is large but not exceptionally large. In half the computer time your calculation needs (for the MCSCF step from what you say) I can do three CASSCF second-order geometry optimization steps followed by a full vibrational calculation and analysis for a calculation with 2.3 million configurations! As I said at the outset, if you post your O2 calculation files (uploaded, not inline please!) we may be able to provide more assistance.

Best regards
Pete

[ankit7540]

Hello Dr. Peter,
Thank you for your comments.
The full output file for O2 was uploaded with the original post. It is given here again.

[taylor]

The fact that there was an output file attached did not come up for me when I looked at the earlier posting, so I was not aware you'd posted it. Sorry about that because I could have replied sooner if I had seen it before!

The statistics presented after the RASSCF calculation converged explain the timing, and this section also (together with the final orbitals and occupation numbers) directly identifies the problem. The calculation is performing vast numbers of microiterations within each macroiteration, and much of this is anyway unproductive because most macroiterations end with some sort of linear dependence among the trial vectors. Convergence of the gradient norm is very slow, although we can see from the statistics that the energy hardly changes in most of the last macroiterations. This explains the timing.

The source of the problem is almost certainly your choice of active space. If you look at the occupation numbers (that they are RAS "pseudo" natural orbitals does not matter in this context) you will see that many of them are small (defined for our purposes here as less than 0.05) and some are very much smaller than this. A crude approximation to the MCSCF gradient and Hessian is given by differences between occupation numbers of the different orbital spaces. This means that the gradient with respect to rotations of (here) the RAS3 orbitals and the secondary orbitals is small, and also that the Hessian matrix elements, which give the curvature of the MCSCF energy with respect to the wave function parameters, will be small. The program is then trying to minimize an energy which is very "flat" with respect to varying the wave function parameters, and convergence will inevitably be very slow. As you see. Another pointer that there is something very unsatisfactory about this wave function choice is that the virial ratio -V/T (printed in each iteration) is about 2.3, when it should be very close to 2, say to within at least two decimal places.

One possibility (although I am doubtful that this is the case) is to start with a more physically plausible set of MOs. I would suggest doing HF before MCSCF, and simply doing a closed-shell calculation, occupation 3 1 1 0 2 0 0 0, so O₂²⁺. The nuclear attraction should produce MOs that are strongly localized in the valence region and which will be a much better starting guess. As I say, this may work, but I am not optimistic.

More realistically, I think you will have to reconsider your active space. How did you choose this in the first place? I notice for instance that the smallest occupation number corresponds to an orbital of delta_u symmetry comprised almost entirely of f character, is this really what you had in mind? I would suggest that rather than trying an extended RASSCF space at the start, you instead construct a valence-shell CASSCF space (I would treat the 2s-derived sigma_g and sigma_u as inactive) and then run a CI calculation, obtaining natural orbitals, and looking at the occupation numbers to see what makes sense to include in a later RASSCF calculation. Orbitals with an occupation number less than about 0.01 are likely to give trouble in such a calculation, and I would probably look to retain only orbitals with an occupation number of 0.03 or larger, at least initially.

Best regards
Pete

[hjaaj]

Hi, I have analyzed your case and identified some problems. And I will provide some comments, and supplement what Pete has written.

note 1) In general I recommend to use MP2 guess for initial orbitals, they are much better than SCF orbitals. My recommended input for your O_2 MCSCF calculation would be (I have marked the lines changed compared to your input with "! changed":

Code: Select all

**WAVE FUNCTIONS
.HF   ! changed
.MP2 ! changed - see note 1
.MCSCF
*SCF INPUT
.DOUBLY
 3 1 1 0 2 0 0 0
.SINGLY
 0 0 0 0 0 1 1 0
*CONFIGURATION INPUT
.SYMMETRY
4
.SPIN MULTIPLICITY
3
.RAS1 SPACE
0 0 0 0 0 0 0 0
.RAS2 SPACE
2 1 1 0 2 1 1 0
.RAS3 SPACE
4 2 2 1 3 2 2 1 ! changed -- see note 2)
.INACTIVE ORBITALS
1 0 0 0 1 0 0 0
.ELECTRONS
12
.RAS1 HOLES
0 0
.RAS3 ELECTRONS
0 2  ! changed -- see note 3)

note 2) You have selected a bad RAS3 active space. Two problems, which are revealed if you look at the natural orbital occupations predicted by MP2: Two very weakly occupied orbitals are included in sym 4 and 8; only one component of the degenerate delta orbitals are included in RAS3. This means that the calculation will break linear symmetry, which is the cause of some of the initial convergence problems. Instead of a RAS3 space of " 4 2 2 2 2 2 2 2" you should have " 4 2 2 1 3 2 2 1".

note 3) with CI-SDT between RAS2 and RAS3 you get the extremely poor convergence of the RASSCF macro iterations which you have observed. With CI-SD you get fast convergence in a few minutes! The reason is basically that with CI-SDT your RASSCF wave function is so close to the CASSCF wave function that the orbital rotations between RAS2 and RAS3 are almost redundant. From mathematics we know that redundancy in the parameters will give very slow convergence - and near redundancy will give slow convergence, as observed for the CI-SDT.

And finally, the integral threshold you modified is not relevant at all for the number of macroiterations, which is what caused the long calculation time.

[ankit7540]

Dear Prof. Taylor and Prof. Jensen,
Thank you for the comments.

I now understand that the RAS3 space chosen is incorrect. Now running with the suggested RAS3 space.

1. The correlation table for D(inf)h to D2h gives the possible symmetries in D2h and in my case I did not chose RAS3 so that the pairs are complete (from delta symmetry). This is, as far as I understand, mentioned in above replies.

2. The reason for choosing a big RAS3 is to have a large virtual space for accounting for polarizability/hyperpolarizability. I am trying to have as good wavefunction possible for accurate response calculations.

3. Process of calculations I have followed is this.

• CASSCF with 21102110 as active space, which finishes fast. (keep the SIRIUS.RST file for next step)

• Multiple RASSCF calculations while incrementally increasing size of RAS3. (Use SIRIUS.RST file for initial orbital data from previous smaller sized calculation)

One more question.

1. My initial guess is that the number of electrons in RAS3 would have effect on polarizability. Its the reason I tried up to 3 e- excitations and also 4e- (not shown here). Are they meaningful ?

Thank you again.

[hjaaj]

re 2 and new question: note that in mcscf the polarizability is to a high degree described with the orbital rotation parameters - a CI-SD as MCSCF will give close to the same accuracy as a CI-SDT with frozen orbitals! When using MCSCF the CI space should be designed to describe the reference state well, properties as polarizability and hyperpolarizabilites are described with the orbital rotation parameters.'

I advice against increasing the RAS3 space beyond what I suggested (this is based on the predicted orbital occupations by MP2). MCSCF is very bad at optimizing orbitals with very low occupation because we use the variational principle and the energy changes very little when replacing an orbital with occupation of e.g. 1.e-4 with one of 2.e-4 (in a full CI). This leads to very slow convergence, because of the mathematics.

Another comment on "one more question", because orbital rotations between RAS2 and RAS3 are included in the RASSCF, the 3- and 4- electron excitations are not important unless they are important for the reference state.

[ankit7540]

Thank you for the insight.

[ankit7540]

Dear Prof. Jensen,
I have tried the suggested dal file for calculation and it finishes very fast. I wonder if we can check the orbital energies after the convergence of the MCSCF wavefunction. Similar to the HF orbital energies (Fock eigenvalues) as following ,

Code: Select all

 Hartree-Fock orbital energies, symmetry 1 ( Ag ),    3 occupied SCF orbitals

       -20.73386803    -1.65134823    -0.73637033     0.14380231     0.20398415
         0.57216090     0.66965868     0.97520603     1.56605162     1.85707752
         2.19462346     2.21163058     3.51254924     4.40483694     6.00825601
         6.08254299     6.88997357     7.97189085    12.63008300

 Hartree-Fock orbital energies, symmetry 2 ( B3u),    1 occupied SCF orbitals

        -0.70616814     0.15463101     0.75904842     0.90541419     1.86980476
         1.92053282     2.69015843     4.11458107     5.93958666     6.31214272
         7.18790577

 Hartree-Fock orbital energies, symmetry 3 ( B2u),    1 occupied SCF orbitals

        -0.70616814     0.15463101     0.75904842     0.90541419     1.86980476
         1.92053282     2.69015843     4.11458107     5.93958666     6.31214272
         7.18790577

 Hartree-Fock orbital energies, symmetry 4 ( B1g),    0 occupied SCF orbitals

         0.57216090     1.85707752     2.21163058     6.08254299     6.88997357

 Hartree-Fock orbital energies, symmetry 5 ( B1u),    2 occupied SCF orbitals

       -20.73292086    -1.10046886     0.10003043     0.32897533     0.50804199
         0.79984785     0.86007323     0.99696502     2.18800661     2.34889843
         2.66967739     2.83863622     2.89803499     6.21084395     6.82044902
         6.93872265     7.04650052     8.29478772    14.36671112

 Hartree-Fock orbital energies, symmetry 6 ( B2g),    0 occupied SCF orbitals

        -0.21187631     0.20268343     0.81804905     1.15816764     2.25787068
         2.32010336     3.10879929     4.94347509     6.50448385     7.46496616
         7.77211108

 Hartree-Fock orbital energies, symmetry 7 ( B3g),    0 occupied SCF orbitals

        -0.21187631     0.20268343     0.81804905     1.15816764     2.25787068
         2.32010336     3.10879929     4.94347509     6.50448385     7.46496616
         7.77211108

 Hartree-Fock orbital energies, symmetry 8 ( Au ),    0 occupied SCF orbitals

         0.79984785     2.18800661     2.83863622     6.93872265     7.04650052

I tried the following,
1. Increasing the PRINT levels, up to 10.

2. Added the following,

Code: Select all

* PRINT LEVELS
.CANONI

to the ** WAVE FUNCTIONS part. In the output, I see line as

Code: Select all

 SIRCAN: canonical/natural orbitals and corresponding
         CI coefficents saved on the Sirius restart file.
 

`SIRIUS.RST` is in binary.

Is there some other way to look at the orbital energies after the convergence. The output file is attached.

p_s1_RAS_O2.out

HF->MP2->MCSCF (RASSCF) for O2

(269.18 KiB) Downloaded 88 times

[hjaaj]

Orbital energies are not well-defined for MCSCF (as they are not well-defined for restricted open-shell Hartree-Fock).

For me, the most important information in analyzing MCSCF results is the natural orbital occupation numbers, they help me decide if I have obtained a balanced and reasonable description of the multi-configuration and/or open-shell nature of the molecule.

[taylor]

What information were you hoping to obtain from "orbital energies" (which, as you have been told, are not well defined for MCSCF)?

Best regards
Pete

[ankit7540]

My plan was to make a plot of energies vs the 8 symmetries and look for distinctions for RAS1, RAS2 and RAS3 spaces.

[taylor]

Then I think you do not need orbital energies because you will get this information from the natural orbital occupation numbers. My usual rules of thumb (your mileage may vary... and others might apply different criteria) are:
If the occupation is
1) Very close to 2, like 1.98x or similar, the orbital should be inactive.
2) Close to 2 like 1.95 or around there, RAS1
3) Close to 0, like 0.05 or thereabouts, RAS3
4) Very close to 0, like 0.02 or less, secondary (i.e., virtual)
5) Anything else, RAS2

Criteria (1) and (4) are an attempt to avoid causing convergence problems with very strongly or very weakly occupied active orbitals.

As the late Bjoern Roos used to say, MCSCF methods are not a black box, and require considerable care and input from the user!

Best regards
Pete

[ankit7540]

Hello Prof. Peter,
Thank you for suggestions regarding orbital occupations. Points (1-3) and (5) are clear to me.
For (4), regarding very small occupation numbers I have a question specific to DALTON.

From the manual I do not see specific parameter for defining these Secondary or Virtual orbitals.

[hjaaj]

This is because the secondary=virtual orbitals by definition are all those you have not put into an inactive or active space - that is, they are simply the rest.