Highlevel workflows
These highlevel workflows are multistep singlepoint energy protocols can either be used on their own as a theorylevel in Singlepoint calculations or used as a SP_theory in workflows such as thermochemprotocol, calc_xyzfiles, confsampler_protocol (see Workflow functionality) or as a theory in run_benchmark (see Benchmarking in ASH) . The ORCA_CC_CBS_Theory uses the ORCA quantum chemistry code for all steps of the workflows and gives a final 0 K electronic energy (no ZPVE). Gradients are not available and these can thus not be used in geometry optimizations or dynamics jobs. The MRCC_CC_CBS_Theory uses the MRCC program (not yet available)
ORCA_CC_CBS_Theory
ORCA_CC_CBS_Theory is synonymous with CC_CBS_Theory.
This is an ASH Theory that carries out a multistep singlepoint protocol to give a CCSD(T)/CBS estimated energy. Multiple ORCA calculations for the given geometry are carried out and the SCF and correlation energies extrapolated to the CCSD(T)/CBS limit using either regular CCSD(T) theory or DLPNOCCSD(T) theory. This workflow is flexible and features multiple ways of approaching the complete basis set limit (CBS) or the complete PNO space limit (CPS). Various options affecting the accuracy, efficiency and robustness of the protocol can be chosen. Many basis set families can be chosen that are available for most of the periodic table. Atomic spinorbit coupling can be automatically included if system is an atom.
class ORCA_CC_CBS_Theory:
def __init__(self, elements=None, scfsetting='TightSCF', extrainputkeyword='', extrablocks='',
guessmode='Cmatrix', memory=5000, numcores=1,
cardinals=None, basisfamily=None, Triplesextrapolation=False, SCFextrapolation=True,
alpha=None, beta=None,
stabilityanalysis=False, CVSR=False, CVbasis="W1mtsmall", F12=False, Openshellreference=None,
DFTreference=None, DFT_RI=False, auxbasis="autoauxmax",
DLPNO=False, pnosetting='extrapolation', pnoextrapolation=[1e6,3.33e7,2.38,'NormalPNO'],
FullLMP2Guess=False, OOCC=False,
T1=False, T1correction=False, T1corrbasis_size='Small', T1corrpnosetting='NormalPNOreduced',
relativity=None, orcadir=None, FCI=False, atomicSOcorrection=False):
Keyword 
Type 
Default value 
Details 


list of strings. 
None 
Required: List of all elements of the molecular system (reaction).
Needed to set up basis set information. Duplicates are OK.
fragment.elems is a valid list.


list of integers 
None 
Required: List of cardinal numbers for basisset extrapolation.
Options: [2,3], [3,4], [4,5] or [5,6]. Singleitem lists also valid:
e.g. [4] (for a single QZ level calculation).


string 
None 
Required: Name of basisset family to use. Various options. See table below. 

string 
None 
Scalar relativity treatment. Options: 'DKH', 'ZORA', 'NoRel', None. 

string 
None 
Path to ORCA. Optional. 

Boolean 
False 
Perform SCF stability analysis for each SCF calculation. 

integer 
1 
Number of cores to use in ORCA calculation. 

Boolean 
False 
Whether to use orbitaloptimized (OOCCD(T)) instead of CCSD(T).
Incompatible with F12 and DLPNO.


string 
None 
Use alternative reference WF in openshell calculation.
Options: 'UHF', 'QRO', None.


string 
None 
Use DFT reference WF (orbitals) in all CC calculations.
Options: (any valid ORCA DFT keyword). Default: None


Boolean 
False 
If using DFTreference, if DFT_RI is True then RIJ/RIJCOSX with
SARC/J and defgrid3 is used to calculate DFT orbitals.


string 
"autoauxmax" 
Auxiliary basis set for CC integrals (/C type). Options: 'autoaux,
'autoauxmax' Default: "autoauxmax"


integer 
5000 
Memory in MBs to use by ORCA. 

string 
'TightSCF' 
SCFconvergence setting in ORCA. Options: 'NormalSCF', 'TightSCF',
'VeryTightSCF', 'ExtremeSCF'.


Boolean 
False 
Use of DLPNO approximation for coupled cluster calculations or not. 

Boolean 
False 
Option to use iterative triples, i.e. DLPNOCCSD(T1) instead of
the default DLPNOCCSD(T0) in all steps.


Boolean 
False 
Option to calculate T1 as a singlestep correction instead. 

string 
'Large' 
Size of basis set in T1 correction. Options: 'Large' (larger cardinal basis),
'Small' (smaller cardinal basis).


string 
'NormalPNOreduced' 
PNO setting for the T1 correction. Options: 'LoosePNO', 'NormalPNO',
'NormalPNOreduced' (TCutPNO=1e6), 'TightPNO'.


Boolean 
False 
Whether to do separate cheaper triples energies extrapolation with
smaller basis sets than singlesdoubles. Requires setting cardinals
to 3 values, e.g. [2,3,4]


list 
[1e6,1e7,1.5,'TightPNO'] 
Parameters for PNOextrapolation (X,Y,Z): X and Y being
TCutPNO thresholds while Z signifies the PNOsetting for the other thresholds.


Boolean 
None 
Whether to use Fulllocal MP2 guess in DLPNO calculations.
Only use if all systems are closedshell.


float 
False 
Manual alpha extrapolation parameter for SCFenergy extrapolation. 

float 
None 
Manual beta extrapolation parameter for correlationenergy extrapolation. 

string 
None 
Optional extra simpleinputkeyword to add in ORCA inputfile. 

string 
None 
Optional extra ORCA blockinput lines to add to ORCA inputfile. 

string 
'CMatrix' 
What ORCA Guessmode to use when doing basisset projections of
orbitals. Options: 'CMatrix' (more robust), 'FMatrix' (cheaper).


Boolean 
False 
Whether to add the experimental atomic spinorbit energy to system
if the system is an atom.


Boolean 
False 
Whether to extrapolate the CCSD(T) calculation to the FullCI limit
by the Goodson formula.


Boolean 
False 
Whether to do explicitly correlated CCSD(T)F12 instead of CCSD(T)/CBS
extrapolation. Use with basisfamily='ccf12'.


Boolean 
False 
Perform additional corevalence+scalarrelativistic correction. 

string 
"W1mtsmall" 
The corevalence basis set to use. The default "W1mtsmall" is only available
for elements HAr. Alternative: some other appropriate corevalence basis set.


Boolean 
True 
Whether the SCF energies are extrapolated or not. If False then the
largest SCF energy calculated will be used (e.g. the def2QZVPP
energy in a def2/[3,4] job).

Basisfamily options
Appropriate allelectron or valence+ECP basis sets for each element with basisfamilies such as : cc, augcc, def2, madef2. If instead an allelectron relativistic approch is desired for all elements then basisfamily="ccdk", "def2zora", "def2dkh" and relativity='DKH' or 'ZORA' can be chosen instead.
Note
"def2" (Ahlrichs allelectron basis sets for HKr, valence basis+def2ECP for KRn)
"madef2" (minimally augmented diffuse Ahlrichs basis sets)
"cc" (correlation consistent basis sets, ccpVnZ for light elements and ccpVnZPP (SKMCDHF ECP) for heavy elements (SrXe, HfRn, Ba, Ru, U)). Note: not available for K.
"augcc" (augmented correlation consistent basis sets, ccpVnZ for light elements and augccpVnZPP for heavy elements)
"ccdk" (DKHrecontracted correlation consistent basis sets, ccpVnZDK for light elements and ccpVnZDK for heavy elements)
"augccdk" (DKHrecontracted aug correlation consistent basis sets, augccpVnZDK for light elements and augccpVnZDK for heavy elements)
"def2zora" (ZORArecontracted Ahlrichs basis sets or SARCZORA basis sets for heavy elements)
"madef2zora" (minimally augmented ZORArecontracted Ahlrichs basis sets or SARCZORA basis sets for heavy elements)
"def2dkh" (DKHrecontracted Ahlrichs basis sets or SARCDKH basis sets for heavy elements)
"def2x2c" (Allelectron X2C relativistic basis sets for HRn)
"madef2dkh" (minimally augmented DKHrecontracted Ahlrichs basis sets or SARCDKH basis sets for heavy elements)
"ccCV" (Corevalence correlation consistent basis sets, ccpwCVnZ)
"augccCV" (augmented corevalence correlation consistent basis sets, augccpwCVnZ)
"ccCVdk" (DKHrecontracted corevalence correlation consistent basis sets, ccpwCVnZDK)
"augccCVdk" (augmented DKHrecontracted corevalence correlation consistent basis sets, augccpwCVnZDK)
"ccCV_3dTMcc_L" (Allelectron DKH protocol for 3d TM complexes. ccpwCVnZDK on 3d transition metals, ccpVNZDK on everything else.)
"augccCV_3dTMcc_L" (Augmented allelectron DKH protocol for 3d TM complexes. ccpwCVnZDK on 3d transition metals, augccpVNZDK on everything else.)
"ccf12" (correlation consistent F12 basis sets for CCSD(T)F12 theory.)
Basisfamily 
Basissets 
Cardinals (n) 
ECP or relativity 

def2 
Ahlrichs def2 on all atoms HRn 

def2ECP on RbRn 
madef2 
Minimally augmented diffuse def2 on all atoms HRn 

def2ECP on RbRn 
def2zora 


relativity='ZORA' 
madef2zora 


relativity='ZORA' 
def2dkh 


relativity='DKH' 
madef2dkh 


relativity='DKH' 
def2x2c 


relativity='DKH' ( later: relativity='X2C') 
cc 


SKMCDHFRSC on SrXe, HfRn, Ba,Ra,U 
ccf12 (use with F12=True) 


SKMCDHFRSC on GaKr, InXe, TlRn 
augcc 


SKMCDHFRSC on SrXe, HfRn, Ba,Ra,U 
ccdk 


relativity='DKH' 
augccdk 


relativity='DKH' 
ccCV 


SKMCDHFRSC on SrXe, HfRn, Ba,Ra,U 
augccCV 


SKMCDHFRSC on SrXe, HfRn, Ba,Ra,U 
ccCVdk 


relativity='DKH' 
augccCVdk 


relativity='DKH' 
ccCV_3dTMcc_L 


relativity='DKH' 
augccCV_3dTMcc_L 


relativity='DKH' 
Note
Note: often missing basis sets for K and Ca. Sometimes there are missing basis sets for specific elements and specific cardinals.
ORCA_CC_CBS_Theory Examples
Basic examples
N2=Fragment(xyzfile='n2.xyz')
cc = ORCA_CC_CBS_Theory(elements=["N"], cardinals = [2,3], basisfamily="cc", numcores=1)
Singlepoint(theory=cc, fragment=N2)
The example above defines an N2 fragment (from file n2.xyz) and runs a singlepoint calculation using the defined ORCA_CC_CBS_Theory object. Multiple CCSD(T) calculations are then carried out using the different basis sets specified by the basisfamily and the cardinals. Cardinals=[2,3] and basisfamily="cc" means that the ccpVDZ and ccpVTZ basis sets will be used. Separate basisset extrapolation of SCF and correlation energies is then performed. Appropriate extrapolation parameters for 2point extrapolations with this basis set family are chosen.
ferrocene=Fragment(xyzfile='ferrocene.xyz')
cc = ORCA_CC_CBS_Theory(elements=["Fe", "C", "H"], cardinals = [2,3], basisfamily="def2", numcores=1,
DLPNO=True, pnosetting="NormalPNO", T1=False)
Singlepoint(theory=cc, fragment=ferrocene)
For a larger molecule like ferrocene, regular CCSD(T) is quite an expensive calculation and so here we invoke the DLPNO approximation via DLPNO=True. We use the 'def2' basis family here with cardinals=[2,3] meaning that the def2SVP and def2TZVPP basis sets will be used. The DLPNO approximation error can be controlled via threshold keywords ('LoosePNO', 'NormalPNO', 'TightPNO'), here we choose 'NormalPNO'. We also choose the regular triples approximation (DLPNOCCSD(T0) by setting T1 to False.
ferrocene=Fragment(xyzfile='ferrocene.xyz')
cc = ORCA_CC_CBS_Theory(elements=ferrocene.elems, cardinals = [3,4], basisfamily="ccCV_3dTMcc_L", relativity='DKH', numcores=1,
DLPNO=True, pnosetting="extrapolation", pnoextrapolation=[6,7] T1=True)
Singlepoint(theory=cc, fragment=ferrocene)
Finally we crank up the accuracy even further by choosing cardinals=[3,4], switch to the basisfamily="ccCV_3dTMcc_L and activate the 'DKH' relativistic approximation. This calculation will utilize a mixed metalligands basis set: ccpwCVTZDK/ccpwCVQZDK on Fe and ccpVDZDK/ccpVTZDK on C,H. Instead of using a single DLPNO threshold we here calculate DLPNOCCSD(T) energies using 2 PNO tresholds and extrapolate to the PNOlimit. Finally we set T1 keyword to True which will tell ORCA to do a more accurate iterative triples DLPNOCCSD(T1) approximation.
For additional examples on using ORCA_CC_CBS_Theory on realworld systems and showing real data see: Tutorial: Highlevel CCSD(T)/CBS workflows
Reaction_Highlevel_Analysis
In order to facilitate the analysis of basisset and/or PNO convergence in CCSD(T) calculations for very simple systems, the Reaction_Highlevel_Analysis function can be used. It will read in an ASH reaction object (containing list of ASH fragments and reaction stoichiometry) and calculate the reaction energy with multiple levels of theory and plot the results using Matplotlib. This allows one to easily see how well converged the results are.
CCSD(T) calculations are performed both with def2 (up to QZ level) and cc basis sets (up to 6Z level), explicitly correlated CCSD(T)F12 calculations (up to QZF12) and complete basis set extrapolations are performed. Note that the largebasis ccpV5Z and ccpV6Z calculations can not be carried out for all systems. Set highest_cardinal to a lower number if required.
Warning
The plots require the Matplotlib library to be installed.
To be added: PNOextrapolation options
def Reaction_Highlevel_Analysis(reaction=None, fraglist=None, stoichiometry=None, numcores=1, memory=7000, reactionlabel='Reactionlabel',
nergy_unit='kcal/mol', extrapolation=True, highest_cardinal=6, plot=True
def2_family=True, cc_family=True, aug_cc_family=False, F12_family=True, DLPNO=False):
"""Function to perform highlevel CCSD(T) calculations for a reaction with associated plots.
Performs CCSD(T) with cc and def2 basis sets, CCSD(T)F12 and CCSD(T)/CBS extrapolations
Args:
reaction ([Reaction object], optional): [ASH Reaction boject]. Defaults to None.
numcores (int, optional): [description]. Defaults to 1.
memory (int, optional): [description]. Defaults to 7000.
def2_family (bool, optional): [description]. Defaults to True.
cc_family (bool, optional): [description]. Defaults to True.
F12_family (bool, optional): [description]. Defaults to True.
highest_cardinal (int, optional): [description]. Defaults to 5.
plot (Boolean): whether to plot the results or not (requires Matplotlib). Defaults to True.
"""
Example (Bond Dissociation Energy of N2):
from ash import *
#Define molecular fragments from XYZfiles or other
N2=Fragment(xyzfile='n2.xyz', charge=0, mult=1, label='N2')
N=Fragment(atom='N', charge=0, mult=4, label='N')
#Define reaction
N2_BDE_reaction = Reaction(fragments=[N2, N], stoichiometry=[1,2], label='N2_BDE', unit='eV')
# Call Reaction_Highlevel_Analysis
Reaction_Highlevel_Analysis(reaction=N2_BDE_reaction, numcores=1, memory=7000,
def2_family=True, cc_family=True, F12_family=True,
extrapolation=True, highest_cardinal=5 )
The outputfile will contain the CCSD(T) total energies and reaction energies for each species and basis set level. Additionally energy vs. basiscardinal plots are created for both the total energy for each species and the reaction energy.
Reaction_FCI_Analysis
With modern approximations to FullCI (selected CI, DMRG, Quantum Monte Carlo etc.) it is possible to obtain a nearFullCI total energy or relative energy that can be used to estimate the accuracy of truncated wavefunction methods (e.g. MP2, CCSD, CCSD(T) etc.). Such an analysis is only possible for relatively small molecules and only for small basis sets, however. ORCA features the ICECI algorithm (a selected CI approach) that can be used for this purpose.
In order to facilitate this kind of analysis ASH features the function Reaction_FCI_Analysis that will automatically run multiple ICECI calculations with ORCA (at userselected thresholds) to estimate the FullCI limit for a given basis set and will then run simpler wavefunction methods with the same basis set for comparison using ORCA. This allows one to see how close e.g. CCSD(T) or CASPT2 is to FullCI for a given energy or relative energy at a specific basis set.
def Reaction_FCI_Analysis(reaction=None, basis=None, basisfile=None, basis_per_element=None,
Do_ICE_CI=True,
MBE_FCI=False, pymbedir=None, mbe_thres_inc=1e5, mbe_orbs_choice='ccsd', mbe_ref_orblist=[],
Do_TGen_fixed_series=True, fixed_tvar=1e11, Do_Tau3_series=True, Do_Tau7_series=True, Do_EP_series=True,
tgen_thresholds=None, ice_nmin=1.999, ice_nmax=0,
separate_MP2_nat_initial_orbitals=True,
DoHF=True,DoMP2=True, DoCC=True, DoCC_CCSD=True, DoCC_CCSDT=True, DoCC_MRCC=False, DoCC_CFour=False, DoCAS=False,
active_space_for_each=None,
DoCC_DFTorbs=True, KS_functionals=['BP86','BHLYP'], Do_OOCC=True, Do_OOMP2=True,
maxcorememory=10000, numcores=1, ice_ci_maxiter=30, ice_etol=1e6,
upper_sel_threshold=1.999, lower_sel_threshold=0,
plot=True, y_axis_label='None', yshift=0.3, ylimits=None, padding=0.4):
Example (Vertical ionization energy of H2O):
from ash import *
#Function to calculate a small molecule reaction energy at the nearFullCI limit at a fixed basis set
#with comparison to simpler methods
#QM code: ORCA
#NearFCI method: ICECI
#Basis set: ccpVDZ
#Molecule: H2O
#Property: VIP
numcores = 1
####################################################################################
#Defining reaction: Vertical IP of H2O
h2o_n = Fragment(databasefile="h2o.xyz", charge=0, mult=1)
h2o_o = Fragment(databasefile="h2o.xyz", charge=1, mult=2)
reaction = Reaction(fragments=[h2o_n, h2o_o], stoichiometry=[1,1], label='H2O_IP', unit='eV')
#What Tgen thresholds to calculate in ICECI?
tgen_thresholds=[5e3,1e3,5e4,1e4,5e5,1e5,5e6]
Reaction_FCI_Analysis(reaction=reaction, basis="ccpVDZ",
Do_Tau3_series=True, Do_Tau7_series=True, Do_TGen_fixed_series=False, fixed_tvar=1e11, Do_EP_series=True,
tgen_thresholds=tgen_thresholds, DoHF=True, DoMP2=True, DoCC=True, maxcorememory=10000, numcores=numcores,
plot=True, y_axis_label='IP', yshift=0.3)
Output:
Warning
The plots require the Matplotlib library to be installed.
WFT theory with flexible orbitalinput option
For many singlereference and multireference methods the reference determinant or orbitals are not optimized, resulting in an inputorbital dependence that can be either mild or sometimes severe. For singlereference CC theory a HF determinant is traditionally used with the assumption that the T1 operator will account indirectly for orbitalrelaxation. There are molecules where this has been found to be a problematic choice with either KohnSham, Brueckner or orbitaloptimized CC references being recommended instead. For selected CI methods such as ICECI (see below), standard CI without orbital optimization is typically performed, resulting in a mild orbital dependency. Typically natural orbitals from a cheaper WF are used in such cases. Natural orbitals are the orbitals that diagonalize the 1electron density matrix.
To facilitate the use of different orbitals in WFT calculations using ORCA, ASH features the ORCA_orbital_setup function which given an orbitalinput choice and basis set runs an ORCA calculation and returns the name of the orbital file to be used for another calculation.
def ORCA_orbital_setup(orbitals_option=None, fragment=None, basis=None, basisblock="", extrablock="", extrainput="",
MP2_density=None, MDCI_density=None, memory=10000, numcores=1, charge=None, mult=None, moreadfile=None,
gtol=2.50e04, nmin=1.98, nmax=0.02, CAS_nel=None, CAS_norb=None,
CASCI=True, tgen=1e4, no_moreadfile_in_CAS=False, ciblockline=""):
The orbitals_option keyword can be an MP2method like: 'MP2', 'RIMP2', 'RISCSMP2', 'OORIMP2' which together with the MP2_density keyword (takes options: 'unrelaxed' or 'relaxed') will result in the calculation of natural orbitals from the chosen MP2 Hamiltonian.
Example: MP2 natural orbitals as input for CCSD(T)
from ash import *
frag = Fragment(databasefile="hf.xyz")
basis="ccpVDZ"
#Using ORCA_orbital_setup to calculate MP2 natural orbitals using the relaxed MP2 density
orbfile, natoccs = ORCA_orbital_setup(fragment=frag, basis=basis, orbitals_option='MP2', MP2_density='relaxed')
#Defining an ORCA CCSD(T) noiter calculation using the natural orbitals from the MP2 calculation
ORCA_CCSD_T = ORCATheory(orcasimpleinput=f"! CCSD(T) {basis} tightscf noiter", moreadfile=orbfile)
Singlepoint(theory=ORCA_CCSD_T, fragment=frag)
Alternatively one can use a CCtype method using the MDCI module in ORCA. Options are: 'CCSD', 'QCISD', 'CEPA/1', 'CPF/1' which together with the MDCI_density keyword (takes options: 'linearized', 'unrelaxed' or 'orbopt') will result in the calculation of natural orbitals using the selected method and density.
Example: QCISD natural orbitals as input for CCSD(T)
from ash import *
frag = Fragment(databasefile="hf.xyz")
basis="ccpVDZ"
#Using ORCA_orbital_setup to calculate QCISD unrelaxed natural orbitals
orbfile, natoccs = ORCA_orbital_setup(fragment=frag, basis=basis, orbitals_option='QCISD', MDCI_density='unrelaxed')
#Defining an ORCA CCSD(T) noiter calculation using the natural orbitals from the QCISD calculation
ORCA_CCSD_T = ORCATheory(orcasimpleinput=f"! CCSD(T) {basis} tightscf noiter", moreadfile=orbfile)
Singlepoint(theory=ORCA_CCSD_T, fragment=frag)
There is also an option to perform a multireference calculation using the CASSCF and MRCI modules in ORCA. Options are: 'CASSCF', 'MRCI', 'MRCI+Q', 'MRSORCI', 'MRDDCI1', 'MRDDCI2', 'MRDDCI3', 'MRAQCC', 'MRACPF'. An active space needs to be defined for this option via CAS_nel and CAS_norb keywords and the user should additionally specify whether CASCI is True or not (if False then CASSCF orbital optimization is carried out). Here the user may want to read in some orbitals via the moreadfile keyword.
Example: MRCI+Q natural orbitals as input for CCSD(T)
from ash import *
frag = Fragment(databasefile="hf.xyz")
basis="ccpVDZ"
#Using ORCA_orbital_setup to calculate MRCI+Q (CAS(2,4) reference) natural orbitals
orbfile, natoccs = ORCA_orbital_setup(fragment=frag, basis=basis, orbitals_option='MRCI+Q',
CAS_nel=2, CAS_norb=4, CASCI=False, moreadfile="orbitals_rotated.gbw")
#Defining an ORCA CCSD(T) noiter calculation using the natural orbitals from the MRCI+Q calculation
ORCA_CCSD_T = ORCATheory(orcasimpleinput=f"! CCSD(T) {basis} tightscf noiter", moreadfile=orbfile)
Singlepoint(theory=ORCA_CCSD_T, fragment=frag)
ICECI workflows
ASH contains a few builtin options to facilitate ICECI or CASCI/CASSCF ICEbased workflows. The function make_ICE_theory allows one to conveniently define ICECI ORCA theories for a given basis set and molecule.
#Create ICECI theory
def make_ICE_theory(basis,tgen, tvar, numcores, nel=None, norb=None, nmin_nmax=False, ice_nmin=None,ice_nmax=None,
autoice=False, basis_per_element=None, maxcorememory=10000, maxiter=20, etol=1e6, moreadfile=None,label=""):
#Workflow to do activespace selection with MP2 or CCSD natural orbitals and then an ICECI based on user thresholds
def Auto_ICE_CAS(fragment=None, basis="ccpVDZ", nmin=1.98, nmax=0.02,
initial_orbitals="MP2", moreadfile=None,
numcores=1, charge=None, mult=None, CASCI=True, tgen=1e4, memory=10000):
AutoICE Example:
Simple way of using the AutoICE option in ORCA. Not necessarily much better than the ORCA way but allows workflows.
from ash import *
numcores=8
frag = Fragment(xyzfile="al2h2_mp2geo.xyz", charge=0, mult=1)
basis="ccpVDZ"
nmin_thresh=1.98
nmax_thresh=0.01
#ICEtheory: based on thresholds
ice = make_ICE_theory("ccpVDZ", 1e4, 1e11,numcores, nmin_nmax=True, ice_nmin=1.98, ice_nmax=0.02, autoice=True,
maxcorememory=10000, label=f"ICE")
result_ICE = ash.Singlepoint(fragment=frag, theory=ice)
Manual AutoICE Example:
Note: Allows manual selection of the natural orbitals calculated and readin. Manually selects the active space based on MP2 occupations and determines the active space which is fed into ICECI calculation.
from ash import *
numcores=8
frag = Fragment(xyzfile="al2h2_mp2geo.xyz", charge=0, mult=1)
basis="ccpVDZ"
nmin_thresh=1.98
nmax_thresh=0.01
#Make MP2 natural orbitals
mp2blocks=f"""
%maxcore 11000
%mp2
natorbs true
density unrelaxed
end
"""
natmp2 = ORCATheory(orcasimpleinput=f"! MP2 {basis} autoaux tightscf", orcablocks=mp2blocks, numcores=numcores, label='MP2', save_output_with_label=True)
Singlepoint(theory=natmp2, fragment=frag)
mofile=f"{natmp2.filename}.mp2nat"
#Determine CAS space based on thresholds
mp2nat_occupations=ash.interfaces.interface_ORCA.MP2_natocc_grab(natmp2.filename+'.out')
print("MP2natoccupations:", mp2nat_occupations)
nel,norb=ash.functions.functions_elstructure.select_space_from_occupations(mp2nat_occupations, selection_thresholds=[nmin_thresh,nmax_thresh])
print(f"Selecting CAS({nel},{norb}) based on thresholds: upper_sel_threshold={nmin_thresh} and lower_sel_threshold={nmax_thresh}")
#ICEtheory: Fixed active space
ice = make_ICE_theory("ccpVDZ", 1e4, 1e11,numcores, nel=nel, norb=norb, maxcorememory=10000, moreadfile=mofile, label=f"ICE")
result_ICE = ash.Singlepoint(fragment=frag, theory=ice)
Manual CASICE Example:
Manual selection of the natural orbitals calculated and readin. Manually selects the active space based on MP2 occupations and determines the active space which is fed into a CASCI/CASSCF calculation using the ICECI algorithm instead of the regular FullCI.
from ash import *
numcores=8
#Input
frag = Fragment(xyzfile="al2h2_mp2geo.xyz", charge=0, mult=1)
basis="ccpVDZ"
nmin_thresh=1.98
nmax_thresh=0.01
tgen=1e4
memory=10000
#Make MP2 natural orbitals
mp2blocks=f"""
%maxcore {memory}
%mp2
natorbs true
density unrelaxed
end
"""
natmp2 = ORCATheory(orcasimpleinput=f"! MP2 {basis} autoaux tightscf", orcablocks=mp2blocks, numcores=numcores, label='MP2', save_output_with_label=True)
Singlepoint(theory=natmp2, fragment=frag)
mofile=f"{natmp2.filename}.mp2nat"
#Determine CAS space based on thresholds
mp2nat_occupations=ash.interfaces.interface_ORCA.MP2_natocc_grab(natmp2.filename+'.out')
print("MP2natoccupations:", mp2nat_occupations)
nel,norb=ash.functions.functions_elstructure.select_space_from_occupations(mp2nat_occupations, selection_thresholds=[nmin_thresh,nmax_thresh])
print(f"Selecting CAS({nel},{norb}) based on thresholds: upper_sel_threshold={nmin_thresh} and lower_sel_threshold={nmax_thresh}")
#ICEtheory: Fixed active space
casblocks=f"""
%maxcore {memory}
%casscf
nel {nel}
norb {norb}
cistep ice
ci
tgen {tgen}
end
end
"""
ice_cas_CI = ORCATheory(orcasimpleinput=f"! CASSCF noiter {basis} tightscf", orcablocks=casblocks, moreadfile=mofile, label=f"ICE")
result_ICE = ash.Singlepoint(fragment=frag, theory=ice_cas_CI)
Automatice CASICE Example:
This is an automatic procedure for the above example but uses CCSD instead of MP2 natural orbitals.
from ash import *
from ash.modules.module_highlevel_workflows import Auto_ICE_CAS
numcores=1
#Fragment
frag = Fragment(xyzfile="al2h2_mp2geo.xyz", charge=0, mult=1)
#Settings
basis="ccpVTZ"
nmin=1.98
nmax=0.02
initial_orbitals="CCSD"
#Call function
Auto_ICE_CAS(fragment=frag, basis=basis, nmin=nmin, nmax=nmax, numcores=numcores, CASCI=True, tgen=1e4, memory=10000,
initial_orbitals=initial_orbitals)
Automatic activespace selection
Similar to above but with more options.
def auto_active_space(fragment=None, orcadir=None, basis="def2SVP", scalar_rel=None, charge=None, mult=None,
initial_orbitals='MP2', functional='TPSS', smeartemp=5000, tgen=1e1, selection_thresholds=[1.999,0.001],
numcores=1):
Workflow to guess a good active space for CASSCF calculation based on a 2step procedure: 1. Calculate MP2natural orbitals (alternative Fractional occupation DFT orbitals) 2. ICECI on top of MP2natural orbitals using a large activespace but with small tgen threshold
Example on ozone:
from ash import *
fragstring="""
O 2.219508975 0.000000000 0.605320629
O 1.305999766 0.913250049 0.557466332
O 2.829559171 0.140210894 1.736132689
"""
fragment=Fragment(coordsstring=fragstrin, charge=0, mult=1)
activespace_dictionary = auto_active_space(fragment=fragment, basis="def2TZVP", charge=0, mult=1,
initial_orbitals='MP2', tgen=1.0)
#Returns dictionary with various active_spaces based on thresholds
Output:
ICECI step done
Note: New natural orbitals from ICECI density matrix formed!
Wavefunction size:
Tgen: 1.0
Tvar: 1e07
Orbital space of CAS(18,37) used for ICECI step
Num generator CFGs: 4370
Num CFGS after S+D: 4370
Table of natural occupation numbers
Orbital MP2natorbs ICEnatocc

0 2.0000 2.0000
1 2.0000 2.0000
2 2.0000 2.0000
3 1.9859 1.9898
4 1.9809 1.9869
5 1.9747 1.9836
6 1.9637 1.9791
7 1.9607 1.9787
8 1.9360 1.9665
9 1.9223 1.9631
10 1.9197 1.9603
11 1.8522 1.9371
12 0.1868 0.0779
13 0.0680 0.0349
14 0.0612 0.0318
15 0.0241 0.0122
16 0.0171 0.0093
17 0.0146 0.0081
18 0.0117 0.0076
19 0.0106 0.0067
20 0.0105 0.0064
...
Recommended active spaces based on ICECI natural occupations:
Minimal (1.95,0.05): CAS(2,2)
Medium1 (1.98,0.02): CAS(12,9)
Medium2 (1.985,0.015): CAS(14,10)
Medium3 (1.99,0.01): CAS(18,13)
Medium4 (1.992,0.008): CAS(18,15)
Large (1.995,0.005): CAS(18,19)
Orbital file to use for future calculations: orca.gbw
Note: orbitals are new natural orbitals formed from the ICECI density matrix