MCSCF#
The core functionality of MultiPsi is multiconfigurational methods, with multiconfigurational self-consistent field (MCSCF) as its cornerstone.
CASSCF#
The most commonly used form of MCSCF nowadays is the complete active space self-consistent field (CASSCF), which is a form of MCSCF where the configurations are generated from a full CI expansion within an active space consisting of a few orbitals. Because of the factorial growth of full CI with the number of active orbitals, the active space should usually be limited to less than 20 electrons in 20 orbitals, with the current record being a CAS(22,22) on the triplet state of hydrogenase. However, this large calculation was only possible on a large cluster using the entire distributed memory available.
To define a CASSCF, we simply need to define our OrbSpace as a full CI expansion (see previous chapters) and provide it to our MCSCF module. For example, we can start with a simple Hartree-Fock, here for furan:
import veloxchem as vlx
import multipsi as mtp
furan_xyz="""9
C -0.86213 -0.90784 0.00007
H -1.63433 -1.64264 -0.00003
C 0.50727 -0.90524 0.00007
C 0.92057 0.47886 -0.00003
C -0.22323 1.23186 -0.00003
O -1.35123 0.40376 -0.00013
H 1.17117 -1.74724 0.00017
H 1.93767 0.81866 0.00007
H -0.46573 2.26986 -0.00013
"""
molecule = vlx.Molecule.from_xyz_string(furan_xyz)
basis = vlx.MolecularBasis.read(molecule,"def2-sv(p)")
scfdrv = vlx.ScfRestrictedDriver()
scf_results = scfdrv.compute(molecule, basis)
Self Consistent Field Driver Setup
====================================
Wave Function Model : Spin-Restricted Hartree-Fock
Initial Guess Model : Superposition of Atomic Densities
Convergence Accelerator : Two Level Direct Inversion of Iterative Subspace
Max. Number of Iterations : 50
Max. Number of Error Vectors : 10
Convergence Threshold : 1.0e-06
ERI Screening Threshold : 1.0e-12
Linear Dependence Threshold : 1.0e-06
* Info * Nuclear repulsion energy: 158.8888494994 a.u.
* Info * Overlap matrix computed in 0.00 sec.
* Info * Kinetic energy matrix computed in 0.00 sec.
* Info * Nuclear potential matrix computed in 0.00 sec.
* Info * Orthogonalization matrix computed in 0.00 sec.
* Info * Starting Reduced Basis SCF calculation...
* Info * ...done. SCF energy in reduced basis set: -228.337309732544 a.u. Time: 0.12 sec.
* Info * Overlap matrix computed in 0.00 sec.
* Info * Kinetic energy matrix computed in 0.00 sec.
* Info * Nuclear potential matrix computed in 0.00 sec.
* Info * Orthogonalization matrix computed in 0.00 sec.
Iter. | Hartree-Fock Energy | Energy Change | Gradient Norm | Max. Gradient | Density Change
--------------------------------------------------------------------------------------------
1 -228.430357327940 0.0000000000 0.16297779 0.01240917 0.00000000
2 -228.433229998540 -0.0028726706 0.05564260 0.00546548 0.06351428
3 -228.433553776746 -0.0003237782 0.00841002 0.00075289 0.01658680
4 -228.433566986233 -0.0000132095 0.00299022 0.00026049 0.00498463
5 -228.433568460938 -0.0000014747 0.00114921 0.00013437 0.00167819
6 -228.433568710663 -0.0000002497 0.00030252 0.00003062 0.00054130
7 -228.433568755137 -0.0000000445 0.00007367 0.00000963 0.00033342
8 -228.433568757939 -0.0000000028 0.00002138 0.00000334 0.00008594
9 -228.433568758030 -0.0000000001 0.00001803 0.00000167 0.00001459
10 -228.433568758068 -0.0000000000 0.00000165 0.00000018 0.00000641
11 -228.433568758069 -0.0000000000 0.00000042 0.00000005 0.00000142
*** SCF converged in 11 iterations. Time: 0.77 sec.
Spin-Restricted Hartree-Fock:
-----------------------------
Total Energy : -228.4335687581 a.u.
Electronic Energy : -387.3224182575 a.u.
Nuclear Repulsion Energy : 158.8888494994 a.u.
------------------------------------
Gradient Norm : 0.0000004225 a.u.
Ground State Information
------------------------
Charge of Molecule : 0.0
Multiplicity (2S+1) : 1
Magnetic Quantum Number (M_S) : 0.0
Spin Restricted Orbitals
------------------------
Molecular Orbital No. 14:
--------------------------
Occupation: 2.000 Energy: -0.58794 a.u.
( 1 C 1p+1: -0.25) ( 3 C 1p+1: 0.27) ( 4 C 1p+1: -0.30)
( 5 C 1p+1: 0.18) ( 5 C 1p-1: -0.18) ( 7 H 1s : 0.25)
( 8 H 1s : -0.25)
Molecular Orbital No. 15:
--------------------------
Occupation: 2.000 Energy: -0.56268 a.u.
( 3 C 1p-1: 0.31) ( 4 C 1p+1: -0.20) ( 4 C 1p-1: -0.24)
( 6 O 1p+1: -0.37) ( 6 O 2p+1: -0.27) ( 7 H 1s : -0.20)
( 8 H 1s : -0.20)
Molecular Orbital No. 16:
--------------------------
Occupation: 2.000 Energy: -0.54763 a.u.
( 1 C 1p+1: -0.25) ( 3 C 1p+1: 0.23) ( 3 C 1p-1: 0.20)
( 4 C 1p-1: -0.29) ( 5 C 1p+1: -0.24) ( 6 O 1p+1: 0.32)
( 6 O 2p+1: 0.24)
Molecular Orbital No. 17:
--------------------------
Occupation: 2.000 Energy: -0.39987 a.u.
( 3 C 1p0 : -0.32) ( 3 C 2p0 : -0.24) ( 4 C 1p0 : -0.32)
( 4 C 2p0 : -0.24) ( 6 O 1p0 : 0.34) ( 6 O 2p0 : 0.30)
Molecular Orbital No. 18:
--------------------------
Occupation: 2.000 Energy: -0.32333 a.u.
( 1 C 1p0 : -0.35) ( 1 C 2p0 : -0.29) ( 3 C 1p0 : -0.23)
( 3 C 2p0 : -0.21) ( 4 C 1p0 : 0.23) ( 4 C 2p0 : 0.21)
( 5 C 1p0 : 0.35) ( 5 C 2p0 : 0.29)
Molecular Orbital No. 19:
--------------------------
Occupation: 0.000 Energy: 0.14744 a.u.
( 1 C 1p0 : -0.32) ( 1 C 2p0 : -0.56) ( 3 C 1p0 : 0.18)
( 3 C 2p0 : 0.35) ( 4 C 1p0 : 0.18) ( 4 C 2p0 : 0.35)
( 5 C 1p0 : -0.32) ( 5 C 2p0 : -0.56) ( 6 O 1p0 : 0.27)
( 6 O 2p0 : 0.39)
Molecular Orbital No. 20:
--------------------------
Occupation: 0.000 Energy: 0.20445 a.u.
( 1 C 3s : -0.65) ( 1 C 2p+1: 0.30) ( 1 C 2p-1: 0.23)
( 2 H 2s : 0.98) ( 3 C 3s : -0.97) ( 3 C 2p+1: -0.29)
( 3 C 2p-1: 0.39) ( 4 C 3s : -0.97) ( 4 C 2p+1: -0.46)
( 4 C 2p-1: -0.17) ( 5 C 3s : -0.65) ( 5 C 2p-1: -0.36)
( 7 H 2s : 1.31) ( 8 H 2s : 1.31) ( 9 H 2s : 0.98)
Molecular Orbital No. 21:
--------------------------
Occupation: 0.000 Energy: 0.22624 a.u.
( 1 C 1p0 : -0.22) ( 1 C 2p0 : -0.54) ( 3 C 1p0 : 0.32)
( 3 C 2p0 : 0.83) ( 4 C 1p0 : -0.32) ( 4 C 2p0 : -0.83)
( 5 C 1p0 : 0.22) ( 5 C 2p0 : 0.54)
Molecular Orbital No. 22:
--------------------------
Occupation: 0.000 Energy: 0.23036 a.u.
( 1 C 3s : 0.73) ( 1 C 2p+1: -0.85) ( 1 C 2p-1: -0.79)
( 2 H 2s : -1.89) ( 3 C 3s : 1.11) ( 3 C 2p-1: 0.24)
( 4 C 3s : -1.11) ( 4 C 2p-1: 0.20) ( 5 C 3s : -0.73)
( 5 C 1p-1: -0.16) ( 5 C 2p+1: 0.28) ( 5 C 2p-1: -1.13)
( 7 H 2s : -0.36) ( 8 H 2s : 0.36) ( 9 H 2s : 1.89)
Molecular Orbital No. 23:
--------------------------
Occupation: 0.000 Energy: 0.23249 a.u.
( 1 C 3s : 0.93) ( 1 C 2p+1: -0.62) ( 1 C 2p-1: -0.57)
( 2 H 2s : -1.56) ( 3 C 3s : -0.60) ( 3 C 2p+1: -0.33)
( 3 C 2p-1: 0.60) ( 4 C 3s : -0.60) ( 4 C 2p+1: -0.60)
( 4 C 2p-1: -0.32) ( 5 C 3s : 0.93) ( 5 C 2p+1: -0.20)
( 5 C 2p-1: 0.81) ( 7 H 2s : 1.14) ( 8 H 2s : 1.14)
( 9 H 2s : -1.56)
Ground State Dipole Moment
----------------------------
X : 0.383475 a.u. 0.974697 Debye
Y : -0.114773 a.u. -0.291725 Debye
Z : 0.000066 a.u. 0.000167 Debye
Total : 0.400283 a.u. 1.017417 Debye
In this example, we will choose all \(\pi\) orbitals (so 5 orbitals total with 6 electrons). The most convenient way to select those orbitals is to use the MultiPsi Orbital Viewer, scroll through the list and click “Active” when we find our desired orbital.
orbviewer=mtp.OrbitalViewer()
orbviewer.plot(molecule,basis,scfdrv.mol_orbs)
We find that the orbitals we want are the number 12, 17, 18, 19 and 21. We can save the resulting orbitals in the file “furan-cas.h5” to use as our starting orbitals.
space=mtp.OrbSpace(molecule,"furan-cas.h5")
mcscfdrv=mtp.McscfDriver()
mc_results = mcscfdrv.compute(molecule,basis,space)
Multi-Configurational Self-Consistent Field Driver
====================================================
Active space definition:
------------------------
Number of inactive (occupied) orbitals: 15
Number of active orbitals: 5
Number of virtual orbitals: 58
This is a CASSCF wavefunction: CAS(6,5)
CI expansion:
-------------
Number of determinants: 100
╭────────────────────────────────────╮
│ Driver settings │
╰────────────────────────────────────╯
Number of states : 1
State-specific calculation
- State of interest : 1
Max. iterations : 50
BFGS window : 5
Convergence thresholds:
- Energy change : 1e-08
- Gradient sq. norm : 1e-08
MCSCF Iterations
-------------------
Iter. | Average Energy | E. Change | Grad. Norm | CI Iter. | Trust rad. | Time
---------------------------------------------------------------------------------
1 -228.473029266 0.0e+00 1.7e-03 0 0.40 0:00:00
2 -228.486533870 -1.4e-02 7.8e-04 0 0.40 0:00:00
3 -228.487831306 -1.3e-03 1.2e-04 0 0.48 0:00:00
4 -228.488000819 -1.7e-04 1.5e-05 0 0.48 0:00:00
5 -228.488013713 -1.3e-05 4.3e-07 0 0.58 0:00:00
6 -228.488014060 -3.5e-07 5.3e-08 0 0.58 0:00:00
7 -228.488014076 -1.6e-08 1.5e-09 0 0.58 0:00:00
8 -228.488014077 -5.6e-10 6.4e-11 0 0.58 0:00:00
** Convergence reached in 8 iterations
9 -228.488014077 -3.4e-11 8.8e-12 0 0.69 0:00:00
Final results
-------------
* State 1
- S^2 : -0.00 (multiplicity = 1.0 )
- Energy : -228.48801407684473
- Natural orbitals
1.99451 1.94212 1.90287 0.09707 0.06342
Spin Restricted Orbitals
------------------------
Molecular Orbital No. 14:
--------------------------
Occupation: 2.000 Energy: -0.55816 a.u.
( 3 C 1p-1: -0.32) ( 4 C 1p+1: 0.19) ( 4 C 1p-1: 0.26)
( 6 O 1p+1: 0.35) ( 6 O 2p+1: 0.26) ( 7 H 1s : 0.20)
( 8 H 1s : 0.20)
Molecular Orbital No. 15:
--------------------------
Occupation: 2.000 Energy: -0.54203 a.u.
( 1 C 1p+1: 0.25) ( 3 C 1p+1: -0.23) ( 3 C 1p-1: -0.18)
( 4 C 1p-1: 0.27) ( 5 C 1p+1: 0.24) ( 6 O 1p+1: -0.34)
( 6 O 2p+1: -0.27)
Molecular Orbital No. 16:
--------------------------
Occupation: 1.995 Energy: -0.58222 a.u.
( 6 O 1p0 : 0.59) ( 6 O 2p0 : 0.43)
Molecular Orbital No. 17:
--------------------------
Occupation: 1.942 Energy: -0.41896 a.u.
( 1 C 1p0 : 0.19) ( 3 C 1p0 : 0.34) ( 3 C 2p0 : 0.23)
( 4 C 1p0 : 0.34) ( 4 C 2p0 : 0.23) ( 5 C 1p0 : 0.19)
( 6 O 2p0 : -0.15)
Molecular Orbital No. 18:
--------------------------
Occupation: 1.903 Energy: -0.30937 a.u.
( 1 C 1p0 : 0.36) ( 1 C 2p0 : 0.28) ( 3 C 1p0 : 0.24)
( 3 C 2p0 : 0.20) ( 4 C 1p0 : -0.24) ( 4 C 2p0 : -0.20)
( 5 C 1p0 : -0.36) ( 5 C 2p0 : -0.28)
Molecular Orbital No. 19:
--------------------------
Occupation: 0.097 Energy: 0.18407 a.u.
( 1 C 1p0 : 0.48) ( 1 C 2p0 : 0.30) ( 3 C 1p0 : -0.30)
( 3 C 2p0 : -0.18) ( 4 C 1p0 : -0.30) ( 4 C 2p0 : -0.18)
( 5 C 1p0 : 0.48) ( 5 C 2p0 : 0.30) ( 6 O 1p0 : -0.28)
( 6 O 2p0 : -0.27)
Molecular Orbital No. 20:
--------------------------
Occupation: 0.063 Energy: 0.28638 a.u.
( 1 C 1p0 : -0.36) ( 1 C 2p0 : -0.24) ( 3 C 1p0 : 0.51)
( 3 C 2p0 : 0.39) ( 4 C 1p0 : -0.51) ( 4 C 2p0 : -0.39)
( 5 C 1p0 : 0.36) ( 5 C 2p0 : 0.24)
Molecular Orbital No. 21:
--------------------------
Occupation: 0.000 Energy: 0.20569 a.u.
( 1 C 3s : -0.62) ( 1 C 2p+1: 0.28) ( 1 C 2p-1: 0.21)
( 2 H 2s : 0.93) ( 3 C 3s : -0.99) ( 3 C 2p+1: -0.30)
( 3 C 2p-1: 0.41) ( 4 C 3s : -0.99) ( 4 C 2p+1: -0.48)
( 4 C 2p-1: -0.18) ( 5 C 3s : -0.62) ( 5 C 2p-1: -0.33)
( 7 H 2s : 1.35) ( 8 H 2s : 1.35) ( 9 H 2s : 0.93)
Molecular Orbital No. 22:
--------------------------
Occupation: 0.000 Energy: 0.23274 a.u.
( 1 C 3s : 0.75) ( 1 C 2p+1: -0.85) ( 1 C 2p-1: -0.79)
( 2 H 2s : -1.89) ( 3 C 3s : 1.12) ( 3 C 2p-1: 0.22)
( 4 C 3s : -1.13) ( 4 C 2p-1: 0.19) ( 5 C 3s : -0.74)
( 5 C 1p-1: -0.16) ( 5 C 2p+1: 0.28) ( 5 C 2p-1: -1.12)
( 7 H 2s : -0.40) ( 8 H 2s : 0.41) ( 9 H 2s : 1.87)
Dipole moment for state 1
---------------------------
X : 0.435398 a.u. 1.106672 Debye
Y : -0.130480 a.u. -0.331646 Debye
Z : 0.000065 a.u. 0.000165 Debye
Total : 0.454529 a.u. 1.155297 Debye
Total MCSCF time: 00:00:02
If the orbitals had changed significantly, we would have had a warning. Still it is good to visualize the final orbitals to check the active space is what we wanted.
orbviewer=mtp.OrbitalViewer()
orbviewer.plot(molecule,basis,space)
It is always recommended to inspect the orbitals after your calculations, as the optimisation may sometimes swap our active orbitals, especially if their occupations are close to 2 or 0.
State-specific and state-averaged CASSCF#
When one is interested in more than just the ground state, for example in spectroscopy or photochemistry, we can specify this in the mcscf driver:
nstates = 3
mc_results = mcscfdrv.compute(molecule,basis,space, nstates)
Multi-Configurational Self-Consistent Field Driver
====================================================
Active space definition:
------------------------
Number of inactive (occupied) orbitals: 15
Number of active orbitals: 5
Number of virtual orbitals: 58
This is a CASSCF wavefunction: CAS(6,5)
CI expansion:
-------------
Number of determinants: 100
╭────────────────────────────────────╮
│ Driver settings │
╰────────────────────────────────────╯
Number of states : 3
State-averaged calculation
- Equal-weights
Max. iterations : 50
BFGS window : 5
Convergence thresholds:
- Energy change : 1e-08
- Gradient sq. norm : 1e-08
MCSCF Iterations
-------------------
Iter. | Average Energy | E. Change | Grad. Norm | CI Iter. | Trust rad. | Time
---------------------------------------------------------------------------------
1 -228.278941396 0.0e+00 1.0e-02 0 0.40 0:00:00
2 -228.296679525 -1.8e-02 1.7e-03 0 0.40 0:00:00
3 -228.298696657 -2.0e-03 2.1e-04 0 0.48 0:00:00
4 -228.298795467 -9.9e-05 1.3e-05 0 0.48 0:00:00
5 -228.298804960 -9.5e-06 5.6e-07 0 0.48 0:00:00
6 -228.298805381 -4.2e-07 3.2e-08 0 0.48 0:00:00
7 -228.298805417 -3.6e-08 2.6e-09 0 0.58 0:00:00
8 -228.298805420 -2.7e-09 2.3e-10 0 0.69 0:00:00
** Convergence reached in 8 iterations
9 -228.298805420 -1.6e-10 1.0e-11 0 0.80 0:00:00
Final results
-------------
* State 1
- S^2 : 0.00 (multiplicity = 1.0 )
- Energy : -228.48189335272252
- Natural orbitals
1.99296 1.94507 1.90727 0.09336 0.06134
* State 2
- S^2 : -0.00 (multiplicity = 1.0 )
- Energy : -228.23493697146856
- Natural orbitals
1.99023 1.43277 1.34328 0.74246 0.49127
* State 3
- S^2 : 0.00 (multiplicity = 1.0 )
- Energy : -228.17958593490408
- Natural orbitals
1.97962 1.95162 1.02709 0.97835 0.06332
Spin Restricted Orbitals
------------------------
Molecular Orbital No. 14:
--------------------------
Occupation: 2.000 Energy: -0.56289 a.u.
( 3 C 1p-1: 0.30) ( 4 C 1p+1: -0.21) ( 4 C 1p-1: -0.22)
( 6 O 1p+1: -0.39) ( 6 O 2p+1: -0.29) ( 7 H 1s : -0.20)
( 8 H 1s : -0.20)
Molecular Orbital No. 15:
--------------------------
Occupation: 2.000 Energy: -0.54801 a.u.
( 1 C 1p+1: -0.25) ( 2 H 1s : 0.15) ( 3 C 1p+1: 0.22)
( 3 C 1p-1: 0.22) ( 4 C 1p-1: -0.31) ( 5 C 1p+1: -0.23)
( 6 O 1p+1: 0.29) ( 6 O 2p+1: 0.22) ( 9 H 1s : 0.15)
Molecular Orbital No. 16:
--------------------------
Occupation: 1.975 Energy: -0.60137 a.u.
( 1 C 1p0 : -0.18) ( 5 C 1p0 : -0.18) ( 6 O 1p0 : -0.54)
( 6 O 2p0 : -0.38)
Molecular Orbital No. 17:
--------------------------
Occupation: 1.762 Energy: -0.37259 a.u.
( 3 C 1p0 : 0.36) ( 3 C 2p0 : 0.26) ( 4 C 1p0 : 0.36)
( 4 C 2p0 : 0.26) ( 6 O 1p0 : -0.23) ( 6 O 2p0 : -0.20)
Molecular Orbital No. 18:
--------------------------
Occupation: 1.412 Energy: -0.24794 a.u.
( 1 C 1p0 : 0.37) ( 1 C 2p0 : 0.30) ( 3 C 1p0 : 0.22)
( 3 C 2p0 : 0.18) ( 4 C 1p0 : -0.22) ( 4 C 2p0 : -0.18)
( 5 C 1p0 : -0.37) ( 5 C 2p0 : -0.30)
Molecular Orbital No. 19:
--------------------------
Occupation: 0.632 Energy: 0.06909 a.u.
( 1 C 1p0 : -0.38) ( 1 C 2p0 : -0.48) ( 3 C 1p0 : 0.18)
( 3 C 2p0 : 0.28) ( 4 C 1p0 : 0.18) ( 4 C 2p0 : 0.28)
( 5 C 1p0 : -0.38) ( 5 C 2p0 : -0.48) ( 6 O 1p0 : 0.33)
( 6 O 2p0 : 0.37)
Molecular Orbital No. 20:
--------------------------
Occupation: 0.220 Energy: 0.21208 a.u.
( 1 C 1p0 : -0.28) ( 1 C 2p0 : -0.38) ( 3 C 1p0 : 0.45)
( 3 C 2p0 : 0.58) ( 4 C 1p0 : -0.45) ( 4 C 2p0 : -0.58)
( 5 C 1p0 : 0.28) ( 5 C 2p0 : 0.38)
Molecular Orbital No. 21:
--------------------------
Occupation: 0.000 Energy: 0.20516 a.u.
( 1 C 3s : 0.63) ( 1 C 2p+1: -0.29) ( 1 C 2p-1: -0.23)
( 2 H 2s : -0.95) ( 3 C 3s : 0.98) ( 3 C 2p+1: 0.30)
( 3 C 2p-1: -0.40) ( 4 C 3s : 0.98) ( 4 C 2p+1: 0.47)
( 4 C 2p-1: 0.17) ( 5 C 3s : 0.64) ( 5 C 2p-1: 0.35)
( 7 H 2s : -1.33) ( 8 H 2s : -1.33) ( 9 H 2s : -0.95)
Molecular Orbital No. 22:
--------------------------
Occupation: 0.000 Energy: 0.23176 a.u.
( 1 C 3s : -0.74) ( 1 C 2p+1: 0.85) ( 1 C 2p-1: 0.79)
( 2 H 2s : 1.88) ( 3 C 3s : -1.11) ( 3 C 2p-1: -0.21)
( 4 C 3s : 1.10) ( 4 C 2p-1: -0.19) ( 5 C 3s : 0.75)
( 5 C 1p-1: 0.16) ( 5 C 2p+1: -0.28) ( 5 C 2p-1: 1.13)
( 7 H 2s : 0.39) ( 8 H 2s : -0.38) ( 9 H 2s : -1.89)
Dipole moment for state 1
---------------------------
X : 0.438507 a.u. 1.114573 Debye
Y : -0.131503 a.u. -0.334248 Debye
Z : 0.000080 a.u. 0.000203 Debye
Total : 0.457801 a.u. 1.163613 Debye
Dipole moment for state 2
---------------------------
X : 0.659697 a.u. 1.676782 Debye
Y : -0.196986 a.u. -0.500688 Debye
Z : -0.000004 a.u. -0.000011 Debye
Total : 0.688479 a.u. 1.749939 Debye
Dipole moment for state 3
---------------------------
X : 0.072580 a.u. 0.184480 Debye
Y : -0.020534 a.u. -0.052193 Debye
Z : -0.000052 a.u. -0.000132 Debye
Total : 0.075429 a.u. 0.191721 Debye
Total MCSCF time: 00:00:02
By default, the program then uses a state-averaged optimisation. This means the orbitals are found to minimize the average of all the selected states (here 3). By doing so, we obtain a compromise and all states are treated at a similar accuracy, favouring energy cancellation in the energy differences (i.e. excitation energies). On the other hand, the orbitals are not optimal for any specific state, and here we can see for example that the ground state energy is slightly higher than the one found in the previous section.
Alternatively, one can decide to optimize the orbitals for a specific state, or biasing the average for a specific state. For this, we use the state-average optimizer but with custom weights instead of the default equal weights:
# no weight on the first and third state, full weight on the second
mcscfdrv.compute(molecule,basis,space, [0,1,0])
Multi-Configurational Self-Consistent Field Driver
====================================================
Active space definition:
------------------------
Number of inactive (occupied) orbitals: 15
Number of active orbitals: 5
Number of virtual orbitals: 58
This is a CASSCF wavefunction: CAS(6,5)
CI expansion:
-------------
Number of determinants: 100
╭────────────────────────────────────╮
│ Driver settings │
╰────────────────────────────────────╯
Number of states : 3
State-specific calculation
- State of interest : 2
Max. iterations : 50
BFGS window : 5
Convergence thresholds:
- Energy change : 1e-08
- Gradient sq. norm : 1e-08
MCSCF Iterations
-------------------
Iter. | Average Energy | E. Change | Grad. Norm | CI Iter. | Trust rad. | Time
---------------------------------------------------------------------------------
1 -228.234936971 0.0e+00 6.4e-03 0 0.40 0:00:00
2 -228.239446513 -4.5e-03 7.4e-04 0 0.40 0:00:00
3 -228.240111313 -6.6e-04 6.5e-05 0 0.48 0:00:00
4 -228.240166877 -5.6e-05 7.9e-06 0 0.48 0:00:00
5 -228.240174439 -7.6e-06 4.4e-07 0 0.58 0:00:00
6 -228.240174649 -2.1e-07 3.7e-08 0 0.58 0:00:00
7 -228.240174665 -1.5e-08 7.7e-10 0 0.58 0:00:00
8 -228.240174665 -5.6e-10 5.7e-11 0 0.58 0:00:00
** Convergence reached in 8 iterations
9 -228.240174665 -4.9e-11 8.5e-12 0 0.69 0:00:00
Final results
-------------
* State 1
- S^2 : -0.00 (multiplicity = 1.0 )
- Energy : -228.4791781877263
- Natural orbitals
1.99448 1.94016 1.89395 0.10619 0.06522
* State 2
- S^2 : 0.00 (multiplicity = 1.0 )
- Energy : -228.24017466518683
- Natural orbitals
1.99260 1.43909 1.31020 0.75252 0.50558
* State 3
- S^2 : 0.00 (multiplicity = 1.0 )
- Energy : -228.1619043945033
- Natural orbitals
1.97985 1.95222 1.02838 0.97599 0.06357
Spin Restricted Orbitals
------------------------
Molecular Orbital No. 14:
--------------------------
Occupation: 2.000 Energy: -0.56225 a.u.
( 3 C 1p-1: 0.35) ( 4 C 1p+1: -0.16) ( 4 C 1p-1: -0.31)
( 6 O 1p+1: -0.28) ( 6 O 2p+1: -0.20) ( 7 H 1s : -0.19)
( 8 H 1s : -0.19)
Molecular Orbital No. 15:
--------------------------
Occupation: 2.000 Energy: -0.54175 a.u.
( 1 C 1p+1: 0.23) ( 3 C 1p+1: -0.22) ( 4 C 1p-1: 0.21)
( 5 C 1p+1: 0.23) ( 6 O 1p+1: -0.41) ( 6 O 2p+1: -0.31)
Molecular Orbital No. 16:
--------------------------
Occupation: 1.993 Energy: -0.59244 a.u.
( 6 O 1p0 : 0.56) ( 6 O 2p0 : 0.42)
Molecular Orbital No. 17:
--------------------------
Occupation: 1.439 Energy: -0.31165 a.u.
( 3 C 1p0 : -0.40) ( 3 C 2p0 : -0.31) ( 4 C 1p0 : -0.40)
( 4 C 2p0 : -0.31) ( 6 O 1p0 : 0.15)
Molecular Orbital No. 18:
--------------------------
Occupation: 1.310 Energy: -0.21315 a.u.
( 1 C 1p0 : -0.43) ( 1 C 2p0 : -0.36) ( 5 C 1p0 : 0.43)
( 5 C 2p0 : 0.36)
Molecular Orbital No. 19:
--------------------------
Occupation: 0.753 Energy: 0.03827 a.u.
( 1 C 1p0 : -0.46) ( 1 C 2p0 : -0.43) ( 3 C 2p0 : 0.21)
( 4 C 2p0 : 0.21) ( 5 C 1p0 : -0.46) ( 5 C 2p0 : -0.43)
( 6 O 1p0 : 0.30) ( 6 O 2p0 : 0.37)
Molecular Orbital No. 20:
--------------------------
Occupation: 0.506 Energy: 0.15207 a.u.
( 1 C 1p0 : -0.18) ( 1 C 2p0 : -0.34) ( 3 C 1p0 : 0.48)
( 3 C 2p0 : 0.63) ( 4 C 1p0 : -0.48) ( 4 C 2p0 : -0.63)
( 5 C 1p0 : 0.18) ( 5 C 2p0 : 0.34)
Molecular Orbital No. 21:
--------------------------
Occupation: 0.000 Energy: 0.20391 a.u.
( 1 C 3s : -0.51) ( 1 C 2p+1: 0.21) ( 1 C 2p-1: 0.15)
( 2 H 2s : 0.74) ( 3 C 3s : -1.05) ( 3 C 2p+1: -0.34)
( 3 C 2p-1: 0.47) ( 4 C 3s : -1.05) ( 4 C 2p+1: -0.54)
( 4 C 2p-1: -0.21) ( 5 C 3s : -0.51) ( 5 C 2p-1: -0.24)
( 7 H 2s : 1.46) ( 8 H 2s : 1.46) ( 9 H 2s : 0.74)
Molecular Orbital No. 22:
--------------------------
Occupation: 0.000 Energy: 0.23646 a.u.
( 1 C 3s : 0.84) ( 1 C 2p+1: -0.73) ( 1 C 2p-1: -0.64)
( 2 H 2s : -1.71) ( 3 C 3s : 1.16) ( 3 C 2p+1: 0.28)
( 4 C 3s : -1.17) ( 4 C 2p+1: -0.28) ( 5 C 3s : -0.83)
( 5 C 2p+1: 0.26) ( 5 C 2p-1: -0.93) ( 7 H 2s : -0.89)
( 8 H 2s : 0.89) ( 9 H 2s : 1.69)
Dipole moment for state 1
---------------------------
X : 0.448404 a.u. 1.139729 Debye
Y : -0.134561 a.u. -0.342021 Debye
Z : 0.000086 a.u. 0.000218 Debye
Total : 0.468159 a.u. 1.189942 Debye
Dipole moment for state 2
---------------------------
X : 0.616427 a.u. 1.566801 Debye
Y : -0.183812 a.u. -0.467205 Debye
Z : 0.000026 a.u. 0.000067 Debye
Total : 0.643249 a.u. 1.634976 Debye
Dipole moment for state 3
---------------------------
X : 0.031888 a.u. 0.081052 Debye
Y : -0.008654 a.u. -0.021996 Debye
Z : 0.000002 a.u. 0.000006 Debye
Total : 0.033042 a.u. 0.083983 Debye
Total MCSCF time: 00:00:02
{'converged': True,
'n_iterations': 9,
'energies': array([-228.47917819, -228.24017467, -228.16190439]),
'ci_vectors': <multipsi.multipsilib.CIVectors at 0x1696418f0>,
'orbital_space': <multipsi.orbspace.OrbSpace at 0x1695df770>,
'ci_expansion': <multipsi.ci_expansion.CIExpansion at 0x169648e60>,
'sa-energy': -228.24017466518683,
'weights': array([0., 1., 0.]),
'Fock': array([[-4.11894154e+01, 4.53009538e-11, -1.97279623e-10, ...,
7.99885989e-09, -3.00240290e-09, -8.01026066e-10],
[ 4.53009839e-11, -2.25245879e+01, 4.23274691e-09, ...,
2.26357663e-09, 2.30278809e-09, -2.19482733e-08],
[-1.97279616e-10, 4.23276737e-09, -2.25245148e+01, ...,
-1.18042254e-08, -1.16745650e-08, -5.01873633e-09],
...,
[ 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, ...,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00],
[ 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, ...,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00],
[ 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, ...,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00]]),
'S2': [-1.6653345369377348e-16, 0.0, 2.220446049250313e-16],
'natural_occupations': [array([1.994484 , 1.94015616, 1.89394575, 0.10619407, 0.06522002]),
array([1.9926043 , 1.43908815, 1.31020459, 0.75252213, 0.50558083]),
array([1.97984814, 1.9522163 , 1.02837774, 0.9759862 , 0.06357163])],
'dipole': {0: array([ 4.48404077e-01, -1.34561349e-01, 8.56111504e-05]),
1: array([ 6.16427096e-01, -1.83812420e-01, 2.62584147e-05]),
2: array([ 3.18881664e-02, -8.65407855e-03, 2.48131865e-06])}}
Here we selected 0 weights on all states except the second one, and as expected, that state is now lower in energy than in the state-averaged case.
Computing each state this way (called the state-specific optimization) provides the best description of each state but suffers from 2 major issues. First, it is not always easy to optimize excited states this way, as the order of the states may change depending on the orbitals and it is thus difficult to truly follow a specific state. Additionally, the orthonormality of the states wavefunction is not enforced, and thus the excited states may appear lower in energy than they should be. For simplicity, for excited states, we recommend using state-averaging or linear response described in a subsequent section.
MC-PDFT#
While MCSCF can recover the correlation within the active space, it does not include any correlation outside of it, which is typically most of the dynamical correlation.
However, in MultiPsi, the MCSCF module also gives us access to MC-PDFT, which could be thought of as an extension of DFT to multiconfigurational wavefunctions. It is enough to provide a functional. The current functionals are “translated” variants of the standard one, indicated by the “t”.
space=mtp.OrbSpace(molecule,"furan-cas.h5")
mcscfdrv=mtp.McscfDriver()
mcscfdrv.xcfun = "tB3LYP"
mc_results = mcscfdrv.compute(molecule,basis,space)
Multi-Configurational Self-Consistent Field Driver
====================================================
Active space definition:
------------------------
Number of inactive (occupied) orbitals: 15
Number of active orbitals: 5
Number of virtual orbitals: 58
This is a CASSCF wavefunction: CAS(6,5)
CI expansion:
-------------
Number of determinants: 100
╭────────────────────────────────────╮
│ Driver settings │
╰────────────────────────────────────╯
Functional : tB3LYP
Global hybrid
- MCSCF fraction : 0.2
Grid level : 4
Number of states : 1
State-specific calculation
- State of interest : 1
Max. iterations : 50
BFGS window : 5
Convergence thresholds:
- Energy change : 1e-08
- Gradient sq. norm : 1e-08
MC-PDFT Iterations
-------------------
Iter. | Average Energy | E. Change | Grad. Norm | CI Iter. | Trust rad. | Time
---------------------------------------------------------------------------------
1 -229.765451621 0.0e+00 1.2e-01 6 0.40 0:00:01
2 -229.785311187 -2.0e-02 1.6e-02 8 0.40 0:00:01
3 -229.789753789 -4.4e-03 3.2e-03 3 0.40 0:00:00
4 -229.790423381 -6.7e-04 1.0e-04 4 0.40 0:00:00
5 -229.790450771 -2.7e-05 9.4e-06 4 0.40 0:00:00
6 -229.790453389 -2.6e-06 2.7e-07 4 0.40 0:00:00
7 -229.790453483 -9.4e-08 5.8e-08 4 0.40 0:00:01
8 -229.790453504 -2.1e-08 3.1e-09 4 0.40 0:00:01
9 -229.790453504 -3.7e-10 1.1e-09 2 0.40 0:00:00
** Convergence reached in 9 iterations
10 -229.790453505 -6.4e-10 2.9e-10 10 0.28 0:00:02
Final results
-------------
* State 1
- S^2 : 0.00 (multiplicity = 1.0 )
- Energy : -229.79045350460603
- Natural orbitals
1.99682 1.98693 1.97701 0.02389 0.01536
Spin Restricted Orbitals
------------------------
Molecular Orbital No. 14:
--------------------------
Occupation: 2.000 Energy: -0.39701 a.u.
( 3 C 1p-1: -0.38) ( 4 C 1p-1: 0.37) ( 6 O 1p+1: 0.18)
( 7 H 1s : 0.18) ( 8 H 1s : 0.18)
Molecular Orbital No. 15:
--------------------------
Occupation: 2.000 Energy: -0.37343 a.u.
( 1 C 1p+1: 0.21) ( 3 C 1p+1: -0.20) ( 4 C 1p+1: -0.16)
( 5 C 1p+1: 0.23) ( 6 O 2s : 0.15) ( 6 O 3s : 0.17)
( 6 O 1p+1: -0.48) ( 6 O 2p+1: -0.32)
Molecular Orbital No. 16:
--------------------------
Occupation: 1.997 Energy: -0.40869 a.u.
( 1 C 1p0 : 0.17) ( 5 C 1p0 : 0.17) ( 6 O 1p0 : 0.58)
( 6 O 2p0 : 0.39)
Molecular Orbital No. 17:
--------------------------
Occupation: 1.987 Energy: -0.30607 a.u.
( 1 C 1p0 : 0.18) ( 3 C 1p0 : 0.37) ( 3 C 2p0 : 0.23)
( 4 C 1p0 : 0.37) ( 4 C 2p0 : 0.23) ( 5 C 1p0 : 0.18)
( 6 O 1p0 : -0.15)
Molecular Orbital No. 18:
--------------------------
Occupation: 1.977 Energy: -0.23154 a.u.
( 1 C 1p0 : -0.39) ( 1 C 2p0 : -0.28) ( 3 C 1p0 : -0.23)
( 3 C 2p0 : -0.17) ( 4 C 1p0 : 0.23) ( 4 C 2p0 : 0.17)
( 5 C 1p0 : 0.39) ( 5 C 2p0 : 0.28)
Molecular Orbital No. 19:
--------------------------
Occupation: 0.024 Energy: 0.02677 a.u.
( 1 C 1p0 : -0.48) ( 1 C 2p0 : -0.28) ( 3 C 1p0 : 0.28)
( 3 C 2p0 : 0.16) ( 4 C 1p0 : 0.28) ( 4 C 2p0 : 0.16)
( 5 C 1p0 : -0.48) ( 5 C 2p0 : -0.28) ( 6 O 1p0 : 0.36)
( 6 O 2p0 : 0.25)
Molecular Orbital No. 20:
--------------------------
Occupation: 0.015 Energy: 0.10091 a.u.
( 1 C 1p0 : -0.36) ( 1 C 2p0 : -0.21) ( 3 C 1p0 : 0.53)
( 3 C 2p0 : 0.37) ( 4 C 1p0 : -0.53) ( 4 C 2p0 : -0.37)
( 5 C 1p0 : 0.36) ( 5 C 2p0 : 0.21)
Molecular Orbital No. 21:
--------------------------
Occupation: 0.000 Energy: 0.08306 a.u.
( 1 C 3s : 0.56) ( 1 C 2p+1: -0.19) ( 1 C 2p-1: -0.27)
( 2 H 2s : -0.85) ( 3 C 3s : 0.83) ( 3 C 2p+1: 0.29)
( 3 C 2p-1: -0.34) ( 4 C 3s : 0.83) ( 4 C 1p+1: 0.18)
( 4 C 2p+1: 0.43) ( 5 C 3s : 0.56) ( 5 C 2p-1: 0.33)
( 7 H 2s : -1.21) ( 8 H 2s : -1.21) ( 9 H 2s : -0.85)
Molecular Orbital No. 22:
--------------------------
Occupation: 0.000 Energy: 0.11203 a.u.
( 1 C 2s : -0.20) ( 1 C 3s : -1.03) ( 1 C 2p+1: 0.58)
( 1 C 2p-1: 0.31) ( 2 H 2s : 1.42) ( 3 C 3s : 0.47)
( 3 C 2p+1: 0.21) ( 3 C 2p-1: -0.42) ( 4 C 3s : 0.47)
( 4 C 2p+1: 0.41) ( 4 C 2p-1: 0.23) ( 5 C 2s : -0.20)
( 5 C 3s : -1.03) ( 5 C 2p+1: 0.31) ( 5 C 2p-1: -0.58)
( 6 O 3s : 0.34) ( 7 H 2s : -0.87) ( 8 H 2s : -0.87)
( 9 H 2s : 1.42)
Dipole moment for state 1
---------------------------
X : 0.283230 a.u. 0.719899 Debye
Y : -0.084824 a.u. -0.215602 Debye
Z : 0.000040 a.u. 0.000101 Debye
Total : 0.295659 a.u. 0.751491 Debye
Total MC-PDFT time: 00:00:12
The full list of available functionals is available by running this:
mtp.XC_Functional.list()
tLDA : Translated LDA functional with Slater exchange and VWN RPA correlation
tLDA5 : Translated LDA functional with Slater exchange and VWN5 correlation
tPBE : Translated PBE exchange and correlation
tPBE-D4 : Translated PBE exchange and correlation with D4 dispersion
tBLYP : Translated BLYP exchange and correlation functionals
tBLYP-D4 : Translated BLYP exchange and correlation functionals with D4 dispersion
tBP86 : Translated BP86 exchange and correlation
tPBE0 : Translated 1-parameter hybrid PBE0 exchange and correlation with (1-l^2) scaled correlation
tPBE0-D4 : Translated 1-parameter hybrid PBE0 exchange and correlation with (1-l^2) scaled correlation with D4 dispersion
tB3LYP : Translated 1-parameter hybrid B3LYP exchange and correlation with (1-l^2) scaled correlation
tB3LYP-D4 : Translated 1-parameter hybrid B3LYP exchange and correlation with (1-l^2) scaled correlation with D4 dispersion
tB3LYP* : Translated 1-parameter hybrid B3LYP* exchange and correlation with (1-l^2) scaled correlation
tB3LYP*-D4 : Translated 1-parameter hybrid B3LYP* exchange and correlation with (1-l^2) scaled correlation with D4 dispersion
tLYP : MCSCF + translated LYP correlation functional
HPG20 : MCSCF + HPG20 correlation functional
tsrLDA : Long-range MCSCF + translated short-range LDA exchange-correlation
tsrPBE : Long-range MCSCF + translated short-range PBE exchange-correlation
tsrBLYP : Long-range MCSCF + translated short-range BLYP exchange-correlation