Complete Active Space Self-Consistent Field
User guide for how to run complete active space self-consistent field (CASSCF) through DaltonProject.
Specifying CAS
To a run a CASSCF calculation through the DaltonProject, the first to do is to specify molecule and basis set:
import daltonproject as dp
molecule, basis = dp.mol_reader('pyridine.mol')
To run a CASSCF calculation, an initial guess of orbitals is needed, one of the simplest guesses is Hartree-Fock orbitals. A Hartree-Fock calculation can be run to generate these orbitals:
hf = dp.QCMethod('HF')
prop = dp.Property(energy=True)
hf_result = dp.dalton.compute(molecule, basis, hf, prop)
print('Hartree-Fock energy =', hf_result.energy)
# Hartree-Fock energy = -243.637999873274
print('Electrons =', hf_result.num_electrons)
# Electrons = 42
No special options are needed for the Hartree-Fock calculation. The simplest possible CAS would be to include the HOMO and LUMO orbitals, i.e. two active orbitals. Since pyridine has 42 electrons, pyridine will have 21 doubly occupied orbitals. If the HOMO is to be included in the active space there will then be 20 inactive orbitals left. This choice will correspond to a CAS(2,2) from canonical Hartree-Fock orbitals. The CASSCF calculation is now ready to be setup.
casscf = dp.QCMethod('CASSCF')
casscf.complete_active_space(2, 2, 20)
casscf.input_orbital_coefficients(hf_result.orbital_coefficients)
prop = dp.Property(energy=True)
casscf_result = dp.dalton.compute(molecule, basis, casscf, prop)
print('CASSCF energy =', casscf_result.energy)
# CASSCF energy = -243.642555771886
The specification complete_active_space(2, 2, 20) is of the form
complete_active_space(active electrons, active orbitals,
inactive orbitals).
The line input_orbital_coefficients(hf_output.orbital_coefficients)
specifies the initial orbitals. Here the initial orbitals are taken from the
previously run Hartree-Fock calculation.
The complete Python script for this example is:
import daltonproject as dp
molecule, basis = dp.mol_reader('pyridine.mol')
hf = dp.QCMethod('HF')
prop = dp.Property(energy=True)
hf_result = dp.dalton.compute(molecule, basis, hf, prop)
print('Hartree-Fock energy =', hf_result.energy)
# Hartree-Fock energy = -243.637999873274
print('Electrons =', hf_result.num_electrons)
# Electrons = 42
casscf = dp.QCMethod('CASSCF')
casscf.complete_active_space(2, 2, 20)
casscf.input_orbital_coefficients(hf_result.orbital_coefficients)
prop = dp.Property(energy=True)
casscf_result = dp.dalton.compute(molecule, basis, casscf, prop)
print('CASSCF energy =', casscf_result.energy)
# CASSCF energy = -243.642555771886
Picking CAS based on Natural Orbital Occupations
The Dalton Project provides tools to help select a CAS based on natural orbital occupations. As an example, let us use these tools on MP2 natural orbital occupations.
import daltonproject as dp
molecule, basis = dp.mol_reader('pyridine.mol')
molecule.analyze_symmetry()
mp2 = dp.QCMethod('MP2')
prop = dp.Property(energy=True)
mp2_result = dp.dalton.compute(molecule, basis, mp2, prop)
print(mp2_result.filename)
# b2773ec9e68d368d0c4b6889a5ec6299edc94101
print('MP2 energy =', mp2_result.energy)
# MP2 energy = -243.9889166118
Note also that symmetry has been enabled. Now the natural orbital occupations can be inspected:
import daltonproject as dp
mp2_output = dp.dalton.OutputParser('b2773ec9e68d368d0c4b6889a5ec6299edc94101')
nat_occs = mp2_output.natural_occupations
print(dp.natural_occupation.scan_occupations(nat_occs))
The method scan_occupations(natural occupation numbers) will per
default list the 14 most important strongly occupied natural orbitals and
weakly occupied natural orbitals. For this example, the print will produce the
output below
Strongly occupied natural orbitals Weakly occupied natural orbitals
Symmetry Occupation Change in occ. Diff. to 2 Symmetry Occupation Change in occ.
2 1.9345 0.0000 0.0655 2 0.0698 0.0000
4 1.9367 0.0022 0.0633 4 0.0630 0.0068
2 1.9674 0.0307 0.0326 2 0.0325 0.0304
3 1.9766 0.0093 0.0234 3 0.0218 0.0107
1 1.9782 0.0016 0.0218 1 0.0207 0.0011
3 1.9788 0.0006 0.0212 3 0.0180 0.0027
1 1.9832 0.0044 0.0168 1 0.0175 0.0005
3 1.9835 0.0003 0.0165 1 0.0166 0.0009
1 1.9862 0.0026 0.0138 3 0.0128 0.0038
1 1.9870 0.0008 0.0130 1 0.0126 0.0002
3 1.9886 0.0016 0.0114 1 0.0119 0.0007
1 1.9895 0.0009 0.0105 3 0.0119 0.0000
3 1.9908 0.0013 0.0092 3 0.0118 0.0001
1 1.9924 0.0015 0.0076 1 0.0105 0.0012
The scan of the natural occupation numbers can aid the user in picking a suitable active space. In this example, the minimal active space has been highlighted in yellow. This active space has been picked as the minimal active space since the most significant jump in natural occupations happens to the next orbitals. The CASSCF calculation can now be setup.
import daltonproject as dp
molecule, basis = dp.mol_reader('pyridine.mol')
molecule.analyze_symmetry()
mp2_output = dp.dalton.OutputParser('b2773ec9e68d368d0c4b6889a5ec6299edc94101')
casscf = dp.QCMethod('CASSCF')
prop = dp.Property(energy=True)
nat_occs = mp2_output.natural_occupations
electrons, cas, inactive = dp.natural_occupation.pick_cas_by_thresholds(nat_occs, 1.94, 0.06)
casscf.complete_active_space(electrons, cas, inactive)
casscf.input_orbital_coefficients(mp2_output.orbital_coefficients)
casscf_result = dp.dalton.compute(molecule, basis, casscf, prop)
print('CASSCF energy =', casscf_result.energy)
# CASSCF energy = -243.699847108852
Since we determined the proper active space, based on natural orbital
occupations, the number of active electrons, active orbitals, and inactive
orbitals does not need to be manually specified. Dalton Project can determine
these quantities based on the threshold we pick for the strongly occupied
natural orbitals and weakly occupied natural orbitals. The method
pick_cas_by_thresholds(nat_occs, 1.94, 0.06) is of the form
pick_cas_by_thresholds(natural occupation numbers,
strongly occupied threshold, weakly occupied threshold).
All strongly occupied natural orbitals below the threshold will be included in
the active space together with all weakly occupied natural orbitals above the
threshold. For this specific example electrons=4,
cas=[0, 2, 0, 2] and inactive=[11, 1, 7, 0].