Choosing an active space for CASSCF using Molpro

Question

I am trying to determine the potential energy curves for a CH2Cl radical by varying the C-Cl distance. For that I use the CASSCF method In MOLPRO with the aug-cc-pVTZ basis set, I didn't get smooth curves. The active space that I choose in Cs symmetry is:

occ(11,3)closed ( 5.1)

I tried other active spaces but I didn't solve the problem so please how can I proceed to solve this problem?

Thanks in advance.

Answer 1

As a developer of MOKIT, I would like to recommend the automatic CASSCF computation using the automr utility in MOKIT. We do not need to manually choose/permute orbitals, no need to specify the active space size, no need to start from a minimal basis set, and no need to take care about the input format of Molpro. All these things can be determined automatically.

Firstly, we write a .gjf file

%mem=48GB
%nprocshared=28
#p CASSCF/aug-cc-pVTZ

mokit{CASSCF_prog=Molpro}

0 2
 C   -0.39943387    0.82617525    0.16329935
 H   -0.04276103    1.33057344   -0.71035216
 H   -1.46943387    0.82618843    0.16329935
Cl    0.36724607    1.91039566    2.04124183

This is just Gaussian gjf format and it can be visualized, e.g. using GaussView. CASSCF_prog=Molpro means you want Molpro to do the CASSCF computation.

Secondly, submit the job

automr CH2Cl.gjf >CH2Cl.out 2>&1

It takes less than 3 minutes and CASSCF is converged in 2 cycles. automr will call Gaussian, GAMESS, Molpro programs successively (in addtion, PySCF is needed), this procedure is actually UHF->UNO->localized paired UNO->GVB->CASSCF. You can find the CASSCF active space size is determined as (3,3) in the output file CH2Cl.out. Related energies and unpaired electrons are printed

E(CASCI)  =      -498.47159238 a.u.
E(CASSCF) =      -498.47282188 a.u.
----------------------- Radical index -----------------------
Not spin singlet. Biradical character will not be computed.
Yamaguchi's unpaired electrons  (sum_n n(2-n)      ):  1.438
Head-Gordon's unpaired electrons(sum_n min(n,(2-n))):  1.230
Head-Gordon's unpaired electrons(sum_n (n(2-n))^2  ):  1.096
-------------------------------------------------------------

The CASSCF initial orbitals and converged orbitals are stored in files

CH2Cl_uhf_uno_asrot2gvb6_s.fch
CH2Cl_uhf_gvb6_CASSCF_NO.fch

respectively, which can be visualized using GaussView, Multiwfn, etc. The CAS(3,3) contains the singly occupied orbital and C-Cl bonding as well as anti-bonding orbitals: (natural orbital occupation numbers are shown below the isosurfaces)

Thirdly, in this step the user has two choices. The first choice is reading CASSCF orbitals of the current geometry and performing Molpro CASSCF calculations of other geometries. In the step above we will obtain files

ch2cl_uhf_gvb6_casscf.a
CH2Cl_uhf_gvb6_CASSCF.com
CH2Cl_uhf_gvb6_CASSCF.xml

The converged CASSCF orbitals are stored in the .xml file. And you can modify the file CH2Cl_uhf_gvb6_CASSCF.com to change geometries. automr is not involved here.

The second choice is again to use automr to perform automatic calculations for each geometry. In this case you may need to specify CASSCF(3,3) in the .gjf file when C-Cl bond is near its equilibrium bond length. This is because a not-stretched C-Cl bond has no multiconfigurational character and thus will not be put into the active space. If you do not specify CASSCF(3,3) in such case, automr may recommend CASSCF(1,1).

Of course, you may think CAS(3,3) is a minimal active space and you want a larger one. This can be easily achieved since all bonding and anti-bonding orbitals are stores in CH2Cl_uhf_uno_asrot2gvb6_s.fch, in which you will also find two C-H bonds, C,Cl core orbitals and Cl valence orbitals. All these orbitals are spatially localized and can be identified with no effort. You can simply perform a CASSCF(13e,10o) calculation using Molpro without permuting orbitals, here the active space includes also the bonding and anti-bonding orbitals of two C-H bonds, three Cl 3s/3p orbitals.

Answer 2

First of all, your syntax isn't exactly the way the MOLPRO manual for the CASSCF program (called "Multi") instructs:

OCC,n1,n2,…,n8; 
CLOSED,n1,n2,…,n8

The following example is from the same page of the MOLPRO documentation:

{multi;occ,9,2;closed,4,1; 

If by (11,3), you mean you have 11 active orbitals with the A' irrep and 3 in the A'' irrep, then I don't find it immediately obvious which 11 orbitals you're choosing.

It might be better for you to try without specifying the occ parameter at all, as the manual says:

"In the absence of an OCC card, the information from the most recent MCSCF calculation is used, or, if there is none, those orbitals corresponding to a minimal valence set, i.e., full valence space, are used."

The "full valence active space" is the one that I would recommend for your first attempt anyway. Before trying this, for a very large number of geometries, I recommend trying it with the STO-3G basis set, then the 6-31G basis set, then the aug-cc-pVDZ basis set, before doing it on the basis set that you've been using, which is aug-cc-pVTZ. The calculations with the smaller basis sets will be a lot faster, so they will burn less of your computer resources. If there's an issue with smoothness of the potential energy curves, then you might want to increase the size of your active space, but if there is no issue, then this "full valence active space" that I'm recommending, has solved your problem!