PySCF interface
PySCF is a very powerful opensource quantum chemistry program (or rather library) with the entire outer interface written in Python and everything else in C, including the very powerful libcint integral library.
ASH features a pretty good interface to PYSCF that allows one to conveniently use the various powerful DFT and WFT based features in the program that can be combined with the geometry optimization, surface scans, NEB, numerical frequencies, QM/MM, MD and metadynamics features of ASH. Due to the nature of PySCF as a Python library (with an essentially unlimited number of options) it is difficult to extensively support every PySCF feature and ASH instead takes the approach of writing wrappercode around the most useful features that makes it suitable for ASH workflows, QM/MM etc. This makes it very easy to use PySCF within ASH for the most basic features but the drawback being that every single PySCF method can not be supported. The pySCF interface in ASH is also used as part of the interface to Block2 and Dice for DMRG, SHCI and QMC calculations. If the use of special features inside PySCF are desired, you may have to use PySCF directly (i.e. outside ASH) or contact us about adding the feature in the ASH interface.
List of features:
One can define DFT, TDDFT, coupled cluster, CASSCF calculations fairly easily with support for various options.
Control over special options such as BSDFT, stability analysis, fractional occupationnoncollinear DFT (GHF/GKS), x2C relativistic Hamiltonians.
Maximum overlap deltaSCF (OODFT) calculations for finding excited SCF states.
Coupled cluster with multiple orbital references. Frozen natural orbital option (FNO).
CASSCF starting orbital options: natural orbitals via MP2, CCSD and CCSD(T). Automatic CAS via DMET_CAS and AVAS.
Read and write checkpointfiles, Molden files and Cubefiles
PySCFTheory interface compatible with BlockTheory and DiceTheory interfaces for DMRG and stochastic heatbath CI calculations.
Support for LOSC PySCF plugin and MCPDFT
Dispersion corrections (D3, D4, TS and MBD) via vdwwrapper
Electrostatically embedded QM/MM including pointcharge gradient is enabled
Limitations:
Not all PySCF features are supported by the interface.
postSCF gradient currently not yet available in the interface
PySCFTheory class:
class PySCFTheory:
def __init__(self, printsetting=False, printlevel=2, numcores=1, label="pyscf",
scf_type=None, basis=None, basis_file=None, ecp=None, functional=None, gridlevel=5, symmetry=False,
guess='minao', dm=None, moreadfile=None, write_chkfile_name='pyscf.chk',
noautostart=False, autostart=True,
soscf=False, damping=None, diis_method='DIIS', diis_start_cycle=0, level_shift=None,
fractional_occupation=False, scf_maxiter=50, direct_scf=True, GHF_complex=False, collinear_option='mcol',
NMF=False, NMF_sigma=None, NMF_distribution=None, stability_analysis=False,
BS=False, HSmult=None,atomstoflip=None,
TDDFT=False, tddft_numstates=10, NTO=False, NTO_states=None,
mom=False, mom_occindex=0, mom_virtindex=1, mom_spinmanifold=0,
dispersion=None, densityfit=False, auxbasis=None, sgx=False, magmom=None,
pe=False, potfile='', filename='pyscf', memory=3100, conv_tol=1e8, verbose_setting=4,
CC=False, CCmethod=None, CC_direct=False, frozen_core_setting='Auto', cc_maxcycle=200, cc_diis_space=6,
CC_density=False, cc_conv_tol_normt=1e06, cc_conv_tol=1e07,
MP2=False,MP2_DF=False,MP2_density=False, DFMP2_density_relaxed=False,
CAS=False, CASSCF=False, CASSCF_numstates=1, CASSCF_weights=None, CASSCF_mults=None,
CASSCF_wfnsyms=None, active_space=None, casscf_maxcycle=200,
frozen_virtuals=None, FNO=False, FNO_orbitals='MP2', FNO_thresh=None, x2c=False,
AVAS=False, DMET_CAS=False, CAS_AO_labels=None, APC=False, apc_max_size=(2,2),
cas_nmin=None, cas_nmax=None, losc=False, loscfunctional=None, LOSC_method='postSCF',
loscpath=None, LOSC_window=None,
mcpdft=False, mcpdft_functional=None):
Keyword 
Type 
Default value 
Details 


string 
None 
Printsetting. if True: printing to standard output otherwise disk. 

string 
None 
Optional label. 

integer 
2 
Printlevel 

integer 
1 
Number of CPU cores used by PySCF 

string 
None 
Type of SCFdeterminant. Options: 'RHF','UHF','RKS','UKS'. 

string 
None 
Name of basis set, can be a string (e.g. 'def2SVP', must be valid PySCF basisname) or a dict with elementspecific keys and valuestrings (basisset name). 

string 
None 
Name of ECP, can be a string ('e.g. 'def2SVP', must be valid PySCF ECPname) or a dict with elementspecific keys and valuestrings (ECP name). 

string 
None 
Name of DFT functional (must be valid PySCF functionalname). 

Boolean 
False 
Use of pointgroup symmetry or not. 

string 
'minao' 
SCF guess options: 'minao', 'atom', 'huckel', 'vsap','1e' 

string 
'' 
Name of DFT functional (must be valid PySCF functionalname). 

Boolean 
False 
Secondorder SCF algorithm active or not. 

float 
None 
Value of damping during SCF. 

string 
'DIIS' 
DIIS method option: 'DIIS', 'ADIIS', 'EDIIS' 

integer 
0 
In which SCF iteration to start the DIIS. 

float 
None 
Value of levelshift (Eh) during SCF. 

Boolean 
False 
Whether fractional occupation is active or not. 

Boolean 
False 
Whether fractional occupation is active or not. 

integer 
50 
Max number of SCF iterations. 

Boolean 
True 
Whether direct SCF algorithm (recalculation of integrals in each iteration) is active or not. False is faster for small systems. 

Boolean 
False 
Whether to use densityfitting (RI) for Coulomb integrals. Use with auxbasis keyword. 

string 
None 
Name of auxiliary basis set to use in densityfitting approximation. Example: 'def2universaljfit'. 

Boolean 
False 
Whether to use seminumerical exchange approximation for HFexchange integrals. Note: gradient is not available 

Boolean 
False 
Whether SCF stability_analysis (calculation of orbital Hessian) is active or not. 

string 
None 
Dispersion correction to use. Options: 'D3', 'D4', 'TS', 'MBD'. Requires pyvdw package. 

string 
None 
Name of PySCF checkpointfile to read in as orbital guess. 

string 
None 
Name of the checkpointfile to write after SCF converges. 

Boolean 
False 
If True, then orbitals are not read in from a checkpointfile. 

list 
None 
If scf_type is 'GHF' or 'GKS', choose magnetic moment: list of the initial collinear spins of each atom. 

Boolean 
False 
If scf_type is 'GHF' or 'GKS', whether complex orbitals are used or not. 

string 
'mcol' 
If scf_type is 'GHF' or 'GKS', collinear option: col, ncol, mcol 

Boolean 
False 
If scf_type is 'GHF' or 'GKS', whether complex orbitals are used or not. 

Boolean 
False 
Whether to find brokensymmetry solution by spinflipping. Requires HSmult, and atomstoflip. 

integer 
None 
BS option: Highspin multiplicity to flip spin from. 

list of integers 
None 
What atom indices to flip spin on. Example: atomstoflip=[0,5] for atoms 0 and 5 

Boolean 
False 
Whether to TDDFT on top of SCF solution or not. 

integer 
10 
Number of TDDFT states calculated. 

Boolean 
False 
Whether to use the X2C scalar relativistic Hamiltonian or not. 

Boolean 
False 
Whether to use a complete active space (CAS) or not. See also CASSCF and active_space keywords below. 

Boolean 
False 
For CAS: Whether CASSCF orbital optimization is active. If False, then CASCI. 

list of integers 
None 
Active space definition (electrons in orbitals), e.g. active_space=[3,2] (3 electrons in 2 orbitals). 

integer 
200 
Maximum number of CASSCF iterations. 

Boolean 
False 
Whether multiconfigurational pair density functional theory (MCPDFT) method is active or not. Requires CAS keywords. 

string 
None 
Name of MCPDFT functional. 

Boolean 
False 
Whether to use the AVAS method to find CAS active space. Requires CAS_AO_labels keyword. 

Boolean 
False 
Whether to use the DMET_CAS method to find CAS active space. Requires CAS_AO_labels keyword. 

list of strings 
None 
List of atomorbital label strings to use in AVAS/DMET_CAS selection. Example: ['Fe 3d', 'Fe 4d', 'C 2pz'] 

float 
None 
If selecting active space from MP2 natural orbitals cas_nmin/cas_nmax tresholds determine active space. 

Boolean 
False 
Whether to use polarizable embedding in PySCF via CPPE library. 

string 
'' 
Name of potential file for in PySCF CPPE polarizable embedding 

string 
'pyscf' 
Filename used for PySCF output 

integer 
3100 
Memory (in MB) used by PySCF . 

float 
1e8 
Convergence tolerance in Eh . 

int 
4 
How verbose PySCF output is. 

Boolean 
False 
Whether to do coupledcluster on top of SCF or not. 

string 
None 
Type of CCSDmethod. Options:'CCSD', 'CCSD(T)'. More options will be available. 

Boolean 
False 
Whether to use integraldirect CC or not. 

integer 
20 
Maximum number of CC iterations. 

string 
'Auto' 
How frozen core is handled. The ASHdefault option is 'Auto' which means that frozen core settings are chosen by ASH (mimics ORCAsettings). 

list 
None 
Optionally freeze selected virtual orbitals in CC calculation. 

Boolean 
False 
Do frozen natural orbital coupled cluster using MP2 natural orbitals. 

float 
None 
Optional threshold to choose virtual natural orbitals to be skipped, based on natural occupation (from MP2 occupations). 

Boolean 
False 
Whether to do localized orbital scaling correction or not. 

string 
None 
The functional used (affect parameters chosen) 

string 
None 
LOSC correction postSCF or full SCF. Options: 'postSCF' or 'SCF' 

list of floats. 
None 
LOSC energy window, e.g. [30,10]. 

string 
None 
Path to losc package. 

Boolean 
False 
Whether to enable the maximum overlap method for deltaSCF calculations. 

integer 
1 
Which relative virtual orbital index to move electron from HOMO into. Default is 1 (LUMO); choose 2 for LUMO+1 etc. 

integer 
0 
What spin manifold to do MOMdeltaSCF calculations in. Default is 0 (i.e. alpha) 
Advanced: PySCFTheory methods
The PySCFTheory class includes several methods that can also be called on their own (if you know what you are doing!)
def create_mol(self, qm_elems, current_coords, charge, mult):
def define_basis(self,basis_string_from_file=None):
def create_mf(self):
def determine_frozen_core(self,elems):
def set_numcores(self,numcores):
def cleanup(self):
def print_orbital_en_and_occ(self,mo_energies=None, mo_occ=None):
def write_orbitals_to_Moldenfile(self,mol, mo_coeffs, occupations, mo_energies=None, label="orbs"):
#Write Cube files for orbital, density or MEP
def cubegen_orbital(self, mol, name, coeffs, nx=60,ny=60,nz=60):
def cubegen_density(self, mol, name, dm, nx=60,ny=60,nz=60):
def cubegen_mep(self, mol, name, dm, nx=60,ny=60,nz=60):
def calculate_natural_orbitals(self,mol, mf, method='MP2', CAS_AO_labels=None, elems=None, relaxed=False, numcores=1):
def calculate_CCSD_natorbs(self,ccsd=None, mf=None):
def calculate_CCSD_T_natorbs(self,ccsd=None, mf=None):
def run_population_analysis(self, mf, unrestricted=True, dm=None, type='Mulliken', label=None, verbose=3):
def run_stability_analysis(self):
def stability_analysis_loop(self,mf,mos,maxcyc=10):
def read_chkfile(self,chkfile):
def setup_guess(self):
def calc_losc(self):
def run_SCF(self,mf=None, dm=None, max_cycle=None):
def run_MP2(self,frozen_orbital_indices=None, MP2_DF=None):
def run_MP2_density(self, mp2object, MP2_DF=None, DFMP2_density_relaxed=None):
def run_CC(self,mf, frozen_orbital_indices=None, CCmethod='CCSD(T)', CC_direct=False, mo_coefficients=None):
def run_CC_density(self,ccobject,mf):
def get_dipole_moment(self, dm=None, label=None):
def get_polarizability_tensor(self):
def set_mf_scfconv_options(self):
def set_mf_smearing(self):
def set_dispersion_options(self):
def set_DF_mf_options(self):
def set_DFT_options(self):
def set_printing_option_mf(self):
def set_collinear_option(self):
def set_frozen_core_settings(self, elems):
def set_embedding_options(self, PC=False):
def density_potential_inversion(self, dm, lambda_par=8, method='ZMP', DF=True):
def run(self, current_coords=None, current_MM_coords=None, MMcharges=None, qm_elems=None,
elems=None, Grad=False, PC=False, numcores=None, pe=False, potfile=None, restart=False, label=None,
charge=None, mult=None):
def prepare_run(self, current_coords=None, current_MM_coords=None, MMcharges=None, qm_elems=None,
elems=None, Grad=False, PC=False, numcores=None, pe=False, potfile=None, restart=False, label=None,
charge=None, mult=None):
def actualrun(self, current_coords=None, current_MM_coords=None, MMcharges=None, qm_elems=None,
elems=None, Grad=False, PC=False, numcores=None, pe=False, potfile=None, restart=False, label=None,
charge=None, mult=None,pyscf=None ):
PySCF installation
The PySCF interface is librarybased and requires a PySCF installation inside the Python environment, typically via Pip (pip install pyscf).
Parallelization
The PySCF parallelization is OpenMP threadbased. The numcores keyword is used to specify the number of threads available to PySCF.
Using the interface
Typicall the pySCFTheory theory object is simply used as an inputtheory object
from ash import *
n2_singlet= Fragment(diatomic="N2", bondlength=1.09, charge=0, mult=1)
#Initialization of the PySCFTheory object (restricted HF here)
pyscf_object = PySCFTheory(basis="ccpVDZ", scf_type='RHF')
#Calling Singlepoint function
Singlepoint(theory=pyscf_object, fragment=n2_singlet)
However, in more advanced usage of the interface you can also call individual methods of the PySCFTheory object. This is considered expertterritory and is typically not recommended.
from ash import *
frag = Fragment(diatomic="N2", bondlength=1.09, charge=0, mult=1)
#Initialization of the PySCFTheory object
pyscf_object = PySCFTheory(basis="ccpVDZ", scf_type='RHF')
#Prepare pySCFTheory object for run: This defines the pyscf mol and mf objects internally
#Also sets various options inside mf and mol object previously defined
pyscf_object.prepare_run(elems=frag.elems, current_coords=frag.coords, charge=frag.charge, mult=frag.mult)
#Setup guess for SCF
pyscf_object.setup_guess()
#Run SCF with optional densitymatrix input (dm) and maxcycle input (here 0, i.e. no SCF)
pyscf_object.run_SCF(dm=None, max_cycle=0) #HFSCF
#Print orbitals, population analysis and dipole
pyscf_object.print_orbital_en_and_occ() #HFSCF
pyscf_object.run_population_analysis(pyscf_object.mf)
pyscf_object.get_dipole_moment()
#Run CC using frozen core
fc_indices=pyscf_object.set_frozen_core_settings(frag.elems)
pyscf_object.run_CC(frozen_orbital_indices=fc_indices, CCmethod='CCSD(T)')
Controlling restart and guess
How an SCF calculations begins can be controlled in a different ways. Internally the SCF guess is handled by the setup_guess method which can be called on its own (see above for example). First it is checked whether the PySCFTheory object already contains a density matrix (dm) and if so, then this is used as the guess. Next it is checked whether the moreadfile keyword has been specified (should contain the name of a pySCF checkpointfile, something.chk) and if so, then the orbitals from the checkpointfile are used as the guess. Next it is checked whether AutoStart has been disabled (either noautostart=True, or autostart=False). Autostart is on by default which means that it will try to read a checkpoint file in the directory with the default filename ("pyscf.chk"). If so then a new orbitalguess is used (based on the guess keyword, defaults to 'minao'). Guess options are: ['minao', 'atom', 'huckel', 'vsap','1e'].
#Reading in a density matrix. some_dm should here be a Numpy array
pyscf_obj = PySCFTheory(scf_type="RHF", basis="def2SVP", dm=some_dm)
#Reading in a checkpoint file using moreadfile
pyscf_obj = PySCFTheory(scf_type="RHF", basis="def2SVP", moreadfile="previous.chk")
#Disabling autostart by autostart=False
pyscf_obj = PySCFTheory(scf_type="RHF", basis="def2SVP", autostart=False)
#Changing guess to huckel
pyscf_obj = PySCFTheory(scf_type="RHF", basis="def2SVP", autostart=False, guess="huckel")
The SCFcontrol functionality above can be utilized to do special things such as performing nonselfconsistent calculations using some energy functional (HF or KSDFT) on some other set of orbitals or density matrix. This requires one to i) read in the orbitals (or the density matrix) and ii) turn off SCF iterations. Performing a nonselfconsistent DFT calculation using HF orbitals/density is called HFDFT (or sometimes densitycorrected DFT, DCDFT) in the literature. An example for this is shown below.
Nonselfconsistent calculation using another set of orbitals (here HFDFT)
#Here we do a HFDFT calculation by running first a HF calculation
#and then using the HF density matrix as a guess for the DFT calculation
from ash import *
frag = Fragment(databasefile="h2o.xyz")
#Run HF calculation from scratch
pySCF_HF = PySCFTheory(scf_type="RHF", basis="def2SVP", autostart=False)
Singlepoint(fragment=frag, theory=pySCF_HF)
#Create DFT object and reading in HF density matrix, also setting scf_maxiter=0 to avoid SCF
pyscf_DFT_HF = PySCFTheory(scf_type="RHF", basis="def2SVP", autostart=False, functional="PBE", dm=pySCF_HF.dm, scf_maxiter=0)
Singlepoint(fragment=frag, theory=pyscf_DFT_HF)
Sometimes in unrestricted SCF calculations, one wants to guide the SCF procedure to find a symmetrybroken solution. This is typically performed in the context of brokensymmetry DFT to describe spincoupled antiferromagnetic states. This can be performed in the PySCF interface by specifying BS=True, setting the spin multiplicity of the highspin state (HSmult) and specifying the atom(s) to flip (atomstoflip is a list of atom indices ).
Brokensymmetry solution via spinflipping a spincenter from the highspin solution
#Here we do a HFDFT calculation by running first a HF calculation
#and then using the HF density matrix as a guess for the DFT calculation
from ash import *
#Specify a BSDFT calculation by setting BS=True and HSmult=3 (highspin multiplicity)
pySCF_HF = PySCFTheory(scf_type="RHF", basis="def2SVP", functional='PBE',
autostart=False, BS=True, HSmult=3, atomstoflip=[0])
Singlepoint(fragment=frag, theory=pySCF_HF, charge=0, mult=1)
Controlling basis set and ECP
Typically it is easiest to specify the basis set as in the examples above using a string (e.g. 'def2SVP'). This option will work as long as the basis set is available inside pySCF and you specify the correct name of the basis set (see pySCF documentation).
If one needs more flexibility, such as using different basis sets for different elements, then one can use a dictionary to specify the basis set for each element.
from ash import *
frag= Fragment(diatomic="HF", bondlength=1.09, charge=0, mult=1)
#Initialization of the PySCFTheory object (restricted HF here)
basis_dict={'H':'ccpVDZ', 'F':'def2SVP'}
pyscf_object = PySCFTheory(basis=basis_dict, scf_type='RHF')
#Calling Singlepoint function
Singlepoint(theory=pyscf_object, fragment=frag)
If one wants to use a custom basis set, e.g. something not available inside pySCF, then one can use the basis_file keyword to specify a file containing the basis set. The basisset file needs to contain the basis set for all elements and needs to be in the NWChem basis set format. It is best to go to the Basis Set Exchange and download the basis set in NWChem format.
from ash import *
frag= Fragment(diatomic="HF", bondlength=1.09, charge=0, mult=1)
#pySCF object using a basisset file
pyscf_object = PySCFTheory(scf_type="UHF", basis_file="bla.basis")
#Calling Singlepoint function
Singlepoint(theory=pyscf_object, fragment=frag)
SCF convergence
In case of SCF convergence problems there are a few options available. One involves modifying the initial guess (see above) or reading in orbitals from a previous calculation (see also above).
If that does not work there are a few other options available such as turning on secondorder SCF (SOSCF), using damping, modifying DIIS startcycle, using levelshifting, enabling fractional occupation as well as increasing max iterations.
Shown below are the relevant keywords with their default values:
PySCFTheory(...,soscf=False, damping=None, diis_method='DIIS', diis_start_cycle=0, level_shift=None,
fractional_occupation=False, scf_maxiter=50)
Controlling integral approximation for Coulomb and HF Exchange
Density fitting for Coulomb and HF Exchange integrals is implemented in pySCF, it is not on by default in the interface. For HF and hybridDFT it is also possible to use seminumerical exchange approximation for HF exchange integrals (similar to RIJCOSX in ORCA).
See https://pyscf.org/user/df.html for more details on what is available in pySCF.
#Density fitting for Coulomb integrals only (recommended for nonhybrid DFT)
#Note: Selecting the efficient Coulombonly auxiliary basis set here
PySCFcalc = PySCFTheory(basis="ccpVDZ", scf_type='RKS', functional='BLYP',
densityfit=True, auxbasis='weigend')
#RIJK: i.e. Density fitting for both Coulomb and HF Exchange (applies if HF or hybrid functional).
#Note: Here we let pySCF automatically choose the RIJK auxiliary basis set (which hopefully exists for the basis set)
PySCFcalc = PySCFTheory(basis="ccpVDZ", scf_type='RKS', functional='BLYP',
densityfit=True)
#Density fitting for Coulomb and + seminumerical Exchange for HF Exchange integrals
#Note: Here choosing again the more efficient Coulombonly auxiliary basis set by Weigend
#Warning: no analytical gradient available for this option
PySCFcalc = PySCFTheory(basis="ccpVDZ", scf_type='RKS', functional='BLYP',
densityfit=False, auxbasis='weigend', sgx=True)
Typical Examples
HFSCF example:
from ash import *
n2_singlet= Fragment(diatomic="N2", bondlength=1.09, charge=0, mult=1)
#Minimal PySCFTheory definitino: RHF calculation
PySCFcalc = PySCFTheory(basis="ccpVDZ", scf_type='RHF')
Singlepoint(theory=PySCFcalc, fragment=n2_singlet)
DFTSCF example:
from ash import *
n2_singlet= Fragment(diatomic="N2", bondlength=1.09, charge=0, mult=1)
#Define PySCF theory: RKSPBE0 hybridDFT calculation
PySCFcalc = PySCFTheory(basis="ccpVDZ", scf_type='RKS', functional="PBE0", gridlevel=6,
numcores=2, memory=3000, filename='pyscf', printsetting=False)
Singlepoint(theory=PySCFcalc, fragment=n2_singlet)
Unrestricted CCSD(T) example:
from ash import *
o2_triplet= Fragment(diatomic="O2", bondlength=1.2075, charge=0, mult=3)
#PySCF with UHF SCF and CCSD(T) on top
PySCFcalc = PySCFTheory(basis="ccpVDZ", numcores=2, scf_type="UHF", CC=True,
CCmethod='CCSD(T)', memory=3000, filename='pyscf', printsetting=False)
Singlepoint(theory=PySCFcalc, fragment=o2_triplet)
Natural orbital calculations from various WF methods
Natural orbitals are defined as the orbitals that diagonalize a oneparticle reduced density matrix (1RDM). If the 1RDM is available from a WF calculation (e.g. MP2, CCSD, CASSCF) then it is possible to calculate the natural orbitals, which is a convenient orbital representation of a manyparticle WF. The naturals orbitals can be visualized (from its MO coefficients) and the natural occupations can be analyzed.
The natural orbitals can be calculated from MP2, CCSD, CCSD(T), CASSCF methods (any method in principle as long as the 1RDM is available). Natural orbitals can also be calculated for DMRG and SHCI WFs (see Block and Dice interfaces for more information). For MP2 and CC methods the 1RDM is not automatically available, you have to request it (e.g. by MP2_density=True or CC_density=True), since it takes additional computational effort. If a CC_density or MP2_density calculation is requested then natural orbitals are automatically calculated and written to a Molden file. Otherwise, the natural orbitals can also be written to a Molden file manually as shown below.
Write Molden files of orbitals
To get access to the orbitals or wavefunction of a pySCF calculation it is easiest to write the orbitals to a Molden file. For WFT calculations it is best to calculate the natural orbitals (see above) and then write the MOcoefficients associated with the natural orbitals to a Molden file. This can be accomplished in a few different ways.
Option 1:
If we are running a simple DFTcalculation (or HF) like below, we can simply call the pySCF_write_Moldenfile function that takes the pySCFTheory object as input and writes the orbitals to a Molden file. Note that the pySCFTheory object needs to have been run (i.e. the Singlepoint calculation is necessary) before.
from ash import *
#Fragment
n2_singlet= Fragment(diatomic="N2", bondlength=1.09, charge=0, mult=1)
#pyscf object
PySCFcalc = PySCFTheory(basis="ccpVDZ", scf_type='RKS', functional="PBE0", gridlevel=6,
numcores=2, memory=3000, filename='pyscf', printsetting=False)
#Singlepoint calculation
Singlepoint(theory=PySCFcalc, fragment=n2_singlet)
#Write moldenfile
pySCF_write_Moldenfile(pyscfobject=PySCFcalc, label="orbs")
Option 2:
Another option is to use the write_orbitals_to_Moldenfile method of the PySCFTheory object. Here we specify the internal mol and mf objects of the PySCFTheory object as arguments to the write_orbitals_to_Moldenfile method. This allows some more flexibility as we could in principle change the input MOcoefficients, MOoccupations and MOenergies (simple numpy arrays). This could e.g. be used to write out natural orbitals if we have diagonalized a correlated density matrix and we have access to the natural orbital MOcoefficients and natural orbital occupations.
from ash import *
n2_singlet= Fragment(diatomic="N2", bondlength=1.09, charge=0, mult=1)
PySCFcalc = PySCFTheory(basis="ccpVDZ", scf_type='RKS', functional="PBE0", gridlevel=6,
numcores=2, memory=3000, filename='pyscf', printsetting=False)
Singlepoint(theory=PySCFcalc, fragment=n2_singlet)
#Using the write_orbitals_to_Moldenfile method of the PySCTheory object
PySCFcalc.write_orbitals_to_Moldenfile(PySCFcalc.mol, PySCFcalc.mf.mo_coeff, PySCFcalc.mf.mo_occ, PySCFcalc.mf.mo_energy, label="orbs")
Option 3:
If you have already run a calculation and you have a PySCF checkpoint file available (only for HF or DFT calculations) then you can create a Moldenfile directly from the checkpoint file using the make_molden_file_PySCF_from_chkfile function (note: you need to import it first). Note that the ASH fragment still needs to be loaded and you have to specify the same basis set as used previously.
from ash import *
from ash.interfaces.interface_pyscf import make_molden_file_PySCF_from_chkfile
#Define fragment again (geometry needed)
n2_singlet= Fragment(diatomic="N2", bondlength=1.09, charge=0, mult=1)
make_molden_file_PySCF_from_chkfile(fragment=n2_singlet, basis="ccpVDZ",
chkfile="pyscf.chk",label="orbs2")
Multireference calculations (CASSCF, MCPDFT etc.)
CASSCF calculations are possible in the interface. Calculations are controlled by the CAS keyword (CAS=True or False) and the CASSCF keyword (CASSCF=True or False). If CAS=True but CASSCF=False then a CASCI calculation is performed (only CI, no orbital optimization). If CAS=True and CASSCF=True then a CASSCF calculation is performed (both CI and orbital optimization). The active space is selected by providing a list of n electrons in m orbitals: active_space=[n,m]. Additionally one can solve for multiple states (controlled by CASSCF_numstates keyword) and it is also possible to specify the multiplicities for each state (CASSCF_mults), weights of the states (CASSCF_weights keyword).
#CASSCF calculation for a singlestate
PySCFcalc = PySCFTheory(basis="ccpVDZ", scf_type='RHF',
CAS=True, CASSCF=True, CASSCF_numstates=1, active_space=[6,5], casscf_maxcycle=200)
#Stateaveraged CASSCF calculations for 3roots with equal weights
PySCFcalc = PySCFTheory(basis="ccpVDZ", scf_type='UHF', CAS=True, CASSCF=True,
CASSCF_numstates=3, active_space=[6,5], CASSCF_mults=[1,3,5], CASSCF_weights=[0.33,0.33,0.33])
A regular HFSCFcalculation is currently automatically performed and can currently not be avoided. However, the HForbital guess for the CASSCF calculation can be controlled in a few different ways. The options are: i) reading in a checkpointfile (moreadfile keyword), ii) use the AVAS automatic active space method (AVAS=True), iii) use the DMET_CAS automatic active space method (DMET_CAS=True), iv) use the APC automatic active space method v) use automatic MP2 natural orbitals.
AVAS and DMET_CAS requires one to set CAS_AO_labels keyword which is a list of atomorbital labels (e.g. ['Fe 3d', 'Fe 4d', 'C 2pz']).
MCPDFT calculations are also possible (mcpdft and mcpdft_functional keywords) but has not been tested.
Excited state calculation examples
TDDFT calculations with NTO analysis
from ash import *
cstring="""
O 0.0 0.0 0.0
H 0.0 0.757 0.587
H 0.0 0.757 0.587
"""
frag = Fragment(coordsstring=cstring, charge=0, mult=1)
pyscf = PySCFTheory(scf_type='RKS', basis='631G', functional='b3lyp',
TDDFT=True, tddft_numstates=10, NTO=True, NTO_states=[1,2])
Singlepoint(theory=pyscf, fragment=frag)
The relevant TDDFT output is shown in the main ASH output like below. Also note that additional output will be present in the pySCF outputfile (by default: pyscf.out)
postSCF is True
Now running TDDFT (Num states: 10)

TDDFT RESULTS

TDDFT transition energies (eV): [ 7.81984875 9.9212168 9.95812916 12.38331843 14.75956804 18.1889349
27.77290941 28.15925452 29.1502703 30.1015163 ]
Transition dipoles: [[2.45304512e01 2.68057788e15 6.69547081e16]
[2.01237402e16 1.21055864e14 6.29424552e01]
[ 2.25211670e15 5.66428336e15 2.04238232e14]
[ 5.34022012e16 5.35950517e01 7.01298803e15]
[ 1.12422599e16 1.06732201e+00 2.04454308e15]
[8.19417866e16 2.28946438e14 7.35926479e01]
[ 3.14351405e14 2.32109432e15 7.66443771e16]
[6.61659079e16 1.36249903e15 1.55571253e01]
[3.49120535e01 2.89921400e15 1.48016888e15]
[5.21074382e15 4.48622759e01 1.51123214e14]]
Oscillator strengths (length): [1.15283538e02 9.62964261e02 1.10832607e28 8.71453664e02
4.11929077e01 2.41342606e01 6.76438154e28 1.66969721e02
8.70464866e02 1.48425612e01]
Oscillator strengths (velocity): [4.05425305e02 1.70258256e01 1.76701202e28 1.13797326e01
3.86383743e01 2.15322494e01 1.72972290e28 1.70834214e02
2.70699046e02 1.03027463e01]
NTO analysis for state 1
Now doing NTO analysis for states: [1, 2]
See pySCF outputfile (pyscf.out) for the NTO analysis
Doing NTO for state: 1
Writing
Doing NTO for state: 2
Writing
pyscf.out contains the following NTO output:
State 1: 7.8198 eV NTO largest component 0.9998830310985499
occNTO: 1.000000 (MO #5)
virNTO: 0.999752 (MO #6)
State 2: 9.92117 eV NTO largest component 0.986473412324918
occNTO: 0.999699 (MO #4)
virNTO: 0.999874 (MO #6)
The NTOorbitals can be visualized using the Moldenfiles created: here ntotd1.molden, ntotd2.molden
deltaSCF calculation using Maximum Overlap Method:
PySCF includes the maximum overlap method that can be used to perform orbitaloptimized SCF calculations of excited states (sometimes called deltaSCF approach). You simply specify the SCFtype, functional and basis set as usual and then specify mom=True and optionally mom_virtindex and mom_spinmanifold keywords.
PySCF will first calculated the groundstate SCF with a regular Aufbau electron configuration and will then modify the guess to move an electron from the HOMO to the specified virtual orbital index (default is mom_virtindex=1 which corresponds to the LUMO) of spinmanifold 0 (alpha). If the SCFtype is restricted (RKS/RHF/ROHF/ROKS) then a ROHF/ROKS calculation will be carried out for the excited SCF calculations. If the SCF type is unrestricted (UKS/UHF) then a UKS/UHF calculation will be carried out.
from ash import *
cstring="""
O 0.0 0.0 0.0
H 0.0 0.757 0.587
H 0.0 0.757 0.587
"""
frag = Fragment(coordsstring=cstring, charge=0, mult=1)
pyscf = PySCFTheory(scf_type='RKS', basis='631G', functional='b3lyp', mom=True, mom_virtindex=1, mom_spinmanifold=0)
Singlepoint(theory=pyscf, fragment=frag)
The output will look like this:

DELTASCF RESULTS

Groundstate SCF energy 76.34781084088975 Eh
Excitedstate SCF energy 76.06068587471486 Eh
deltaSCF transition energy 7.812957454584829 eV
Alpha electron occupation pattern of ground state : [1. 1. 1. 1. 1. 0. 0. 0. 0. 0. 0. 0. 0.]
Beta electron occupation pattern of ground state : [1. 1. 1. 1. 1. 0. 0. 0. 0. 0. 0. 0. 0.]
Alpha electron occupation pattern of excited state : [1. 1. 1. 1. 0. 1. 0. 0. 0. 0. 0. 0. 0.]
Beta electron occupation pattern of excited state : [1. 1. 1. 1. 1. 0. 0. 0. 0. 0. 0. 0. 0.]
deltaSCF calculation using Maximum Overlap Method: