For isolated molecules of modest size the ab initio methods have reached great accuracy at present both for ground and excited states. Theoretical studies on isolated molecules, however, may have limited value to bench chemists since most of the actual chemistry takes place in a solvent. If solute-solvent interactions are strong they may have a large impact on the electronic structure of a system and then on its excitation spectrum, reactivity, and properties. For these reasons, numerous models have been developed to deal with solute-solvent interactions in ab initio quantum chemical calculations. A microscopic description of solvation effects can be obtained by a supermolecule approach or by combining statistical mechanical simulation techniques with quantum chemical methods. Such methods, however, demand expensive computations. By contrast, at the phenomenological level, the solvent can be regarded as a dielectric continuum, and there are a number of approaches [47,48] based on the classical reaction field concept.
MOLCAS incorporates a Self-Consistent Reaction Field (SCRF) model within the framework of the SCF and RASSCF programs. The reaction field formalism is based on a molecule placed in a cavity of a dielectric medium. The surrounding is characterized by its electric permitivity and a radius describing a cavity (containing all the atoms of the studied molecule) indicating where the dielectric medium starts. The charge distribution of the molecule will introduce an electric field acting on the dielectric medium. This reaction field will interact with the molecule. The interaction will manifest itself as a perturbation to the one-electron Hamiltonian. The perturbation will be automatically computed in a direct fashion and added to the one-electron Hamiltonian. Once the reaction field has been computed by the SCF or RASSCF programs it can be added as a perturbation to the one-electron Hamiltonian at all levels of theory, such as MRCI or CASPT2.
The present version of the model only uses spherical cavities. In addition, it includes Pauli repulsion due to the medium by introducing a repulsive potential representing the exchange repulsion between the solute and the solvent. This is done by defining a penalty function of gaussian type, generating the corresponding spherical well integrals, and adding them to the one-electron Hamiltonian. When the repulsion potential is used, the size of the cavity should be optimized for the ground state of the molecule (see below). If the repulsive potential is not used and the cavity size is chosen to be smaller (molecular size plus van der Waals radius as is the usual choice in the literature) one must be aware of the consequences: larger solvent effects but also an unknown presence of molecular charge outside the boundaries of the cavity. This is not a consequence of the present model but it is a general feature of cavity models [48].
Dielectric cavity models, in general, assume equilibrium between the electronic state of the solute and the reaction field. This condition is not fulfilled for an electronic excitation. Therefore, to compute excited states MOLCAS introduces the time dependence into the RASSCF program by partitioning the reaction field factor into two parts, a slow and a fast component. The slow field follows geometrical changes of the solute and is, in practice, determined by the properties of the ground state. The fast component can be considered as the instantaneous electronic polarization that follows the absorption of a photon and the strength of this fast component is proportional to the square of the refractive index, that is, the optical value of the dielectric constant [48].
We begin by performing a CASSCF/CASPT2 reaction field calculation on the ground state of a molecule. The reaction field calculations are performed in MOLCAS by introducing the reaction field parameters in the input to the integral program SEWARD via the input keyword
RF-Input
If no repulsive potential is going to be used the input simply consist in adding the above mentioned keyword including the appropriate data (dielectric constant of the medium, cavity size, and angular quantum number of the highest multipole moment of the charge distribution) into the SEWARD input:
&SEWARD &END ... ... RF-Input Reaction field 80 8.0 4 End of RF-Input ... ... End of Input
This will compute the reaction field at those levels. The dielectric constant 80.0 correspond to water as solvent. The radius of the cavity is 8.0 in atomic units. Finally 4 is the maximum angular moment number used in the multipole expansion. The cavity origin is the coordinate origin, thus the molecule must be placed accordingly. If we want to include the reaction field at other levels of theory the keyword RFPErt must be added to the MOTRA or CASPT2 inputs.
We are, however, going to explain the more complicated situation where a repulsive well potential has to be added to the model. In this case it is convenient to optimize the size of the cavity, although in so doing we obtain large cavity sizes and therefore smaller solvent effects. More realistic results can be obtained if additional and specific solvent molecules are added inside the cavity.
To define the well potential we have to add the keyword WELL Integrals to the SEWARD input to compute and add the Pauli repulsion integrals to the bare Hamiltonian.
The requirements considered to build this potential
are that it shall reproduce solvation energies for spherical particles,
ions, and that it must be wide enough so that the electrons in the
excited state of the molecules are also confined to the cavity. Negative ions
have the property that their electrons are loosely bound and they are thus suited
for parametrizing the repulsive potential.
The final result of different calibration calculations
[49,48] is a penalty function which includes four
Gaussians. If 
is
the radius of the cavity the Gaussians are placed at distances 
,
, 
and 
a.u. from the cavity's center with
exponents 5.0, 3.5, 2.0 and 1.4, respectively.
As an example we will use the N,N-dimethylaminobenzonitrile (DMABN) molecule (see Figure 4.6). This is a well known system with large dipole moments both in ground and excited states which suffer important effects due to the polar environment.
&SEWARD &END
Title
para-DMABN molecule. Cavity size: 10 au.
Symmetry
X XY
Basis set
N.ANO-S...3s2p1d.
N1 0.0000000000 0.0000000000 4.7847613288
N2 0.0000000000 0.0000000000 -8.1106617786
End of basis
Basis set
C.ANO-S...3s2p1d.
C1 0.0000000000 0.0000000000 2.1618352923
C2 0.0000000000 2.2430930886 0.7747833630
C3 0.0000000000 2.2317547910 -1.8500321252
C4 0.0000000000 0.0000000000 -3.1917306021
C5 0.0000000000 0.0000000000 -5.9242601761
C6 0.0000000000 2.4377336900 6.0640991723
End of basis
Basis set
H.ANO-S...2s.
H1 0.0000000000 4.0043085530 -2.8534714086
H2 0.0000000000 4.0326542950 1.7215314260
H3 0.0000000000 2.1467175630 8.0879851846
H4 1.5779129980 3.6622699270 5.5104123361
End of basis
RF-Input
reaction field
38.8 10.0 4
End of RF-Input
Well Int
4
1.0 5.0 12.0
1.0 3.5 13.0
1.0 2.0 15.0
1.0 1.4 17.0
End of Input
&SCF &END
TITLE
DMABN molecule
OCCUPIED
20 2 12 5
ITERATIONS
50
END OF INPUT
&RASSCF &END
TITLE
p-DMABN
SYMMETRY
1
SPIN
1
NACTEL
10 0 0
FROZEN
8 0 3 0
INACTIVE
12 1 9 1
RAS2
0 2 0 7
THRS
1.0E-06,1.0E-03,1.0E-03
ITER
50,25
LUMORB
END OF INPUT
In the SEWARD input the WELL Integrals must include
first the number of gaussians used (four), followed by the
coefficient and exponent of the gaussian and the radius of
the cavity in the sequence explained above: first the most
compact gaussian with the radius plus 2.0 au, and so on
to the least compact gaussian.
Here, we have defined a cavity size of 10 au
(cavity centered at coordinate origin). The RASSCF program will read
the RCTFLD input, prepared this time for acetonitrile
(
= 38.8 ), a cavity size of 10.0 au (the same
as in the SEWARD input) and a multipole expansion up to the fourth order
which is considered sufficient [48].
The active space includes the 
space over
the molecular plane, excluding the orbital of the CN
group which lies in the molecular plane.
We repeat the calculation for different cavity sizes in order to find the radius which gives the lowest absolute energy at the CASSCF level. The presence of the repulsive terms allows the cavity radius to be computed by energy minimization. For the calculations using different cavity sizes it is not necessary to repeat the calculation of all the integrals, just those related to the well potential. Therefore, the keyword ONEOnly can be included in the SEWARD input. The ONEINT file will be modified and the ORDINT file is kept the same for each molecular geometry. The energies obtained are in Table 4.13.
|
Taking the gas-phase value (no cav.) as the reference, the CASSCF energy obtained with a 10.0 au cavity radius is higher. This is an effect of the repulsive potential, meaning that the molecule is too close to the boundaries. Therefore we discard this value and use the values from 11.0 to 16.0 to make a simple second order fit and obtain a minimum for the cavity radius at 13.8 au.
Once we have this value we also need to optimize the position of the
molecule in the cavity. Some parts of the molecule, especially those
with more negative charge, tend to move close to the
boundary. Remember than the sphere representing the cavity has
its origin in the cartesian coordinates origin. We use the radius of
13.8 au and compute the CASSCF energy at different displacements
along the coordinate axis. Fortunately enough, this molecule has
C
symmetry. That means that displacements along two of the
axis (
and 
) are restricted by symmetry. Therefore it is
necessary to analyze only the displacements along the 
coordinate.
In a less symmetric molecule all
the displacements should be studied even including combination of the displacements.
The result may even be a three dimensional net, although no
great accuracy is really required. The results for DMABN n 
symmetry are compiled
in Table 4.14.
|
Fitting these values to a curve we obtain an optimal displacement of -1.0 au. We move the molecule and reoptimize the cavity radius at the new position of the molecule. The results are listed in Table 4.15.
|
There is no significant change. The cavity radius is then selected as 13.8 au and the position of the molecule with respect to the cavity is kept as in the last calculation. The calculation is carried out with the new values. The SCF or RASSCF outputs will contain the information about the contributions to the solvation energy. The CASSCF energy obtained will include the reaction field effects and an analysis of the contribution to the solvation energy for each value of the multipole expansion:
Reaction field specifications:
------------------------------
~
Dielectric Constant : .388E+02
Radius of Cavity(au): .138E+02
Truncation after : 4
~
Multipole analysis of the contributions to the dielectric solvation energy
~
--------------------------------------
l dE
--------------------------------------
0 .0000000
1 -.0013597
2 -.0001255
3 -.0000265
4 -.0000013
--------------------------------------
As was explained above, the calculation of the solvent effects on excited states must be done in a different way. The inclusion of the solvent effects on the ground state of one molecule leads to a relaxation both for nuclei and electrons. It is considered to be in equilibrium with the environment. However, when an absorption or emission occurs, a dynamic process starts. It is usually considered that the new electronic distribution of the final state polarizes the solvent, which affects again in a different way the electronic part of the molecule than in the ground state case. This is considered to be a fast motion, where the equilibrium is achieved by solvent and solute electrons. The situation is not the same for nuclei. Slow motions are included in this case, and the nuclei do not achieve the equilibrium with solvent. These slow motions and the related interactions remain as in the ground state. Therefore, the calculations must be carried out using this sequence of steps:
HFrctfld 1 1 1.796
The first number, one, indicates that this is a ground (or initial) state calculation and the slow part or the reaction field is stored in the file COMFILE. The second number is the number of the computed root, one as we are interested in the ground state at this stage. The third number is the square of the refractive index of the solvent, also called the dielectric constant at high frequencies or optical value (from an approximate expression of the Clausius-Mossoti formulae). The RF-Input section in the integral input must be included also, as in the previous case :
&SEWARD &END ... ... RF-Input reaction field 38.8 13.8 4 end of RF-Input ... ... End of Input
-pole
component. The slow part is calculated for the ground state wave function
as the (
/
) fraction of the reaction field,
where is the reaction field factor
 |
(4.5) |
is the reaction field
factor using
, 
is the refraction
index for some suitable frequency in the infrared or visible range,
and 
is the component 
of a -pole moment of order 
.
This value, 
, is stored in the communications
file together with the corresponding
self-energy used to create the cavity. These fields in the COMFILE file
will not be replaced until we perform another RASSCF calculation using the
same COMFILE and the first value in the keyword HFRCtfld
set to one.
HFrctfld 0 2 1.796
where the zero indicates this is an excited state calculation, the two is the number of the computed root, and 1.796 the square of the refractive index of the solvent. The same RCTFLD namelist is used as before.
The total reaction field of the excited state will be added to the COMFILE file, which can be then used as a perturbation on other program for that particular excited state. The slow part of the reaction field of the initial state remains in the COMFILE.
Remember, however, that each time an SCF or RASSCF calculation is done the resulting total reaction field is stored in the communications file (COMFILE). Therefore, one must control which file is being used. The following is an easy sequence for the calculations:
&SEWARD &END
Title
para-DMABN molecule. Cavity size: 13.8 au.
Symmetry
X XY
Basis set
N.ANO-S...3s2p1d.
N1 0.0000000000 0.0000000000 3.7847613288
N2 0.0000000000 0.0000000000 -9.1106617786
End of basis
Basis set
C.ANO-S...3s2p1d.
C1 0.0000000000 0.0000000000 1.1618352923
C2 0.0000000000 2.2430930886 -0.2252166370
C3 0.0000000000 2.2317547910 -2.8500321252
C4 0.0000000000 0.0000000000 -4.1917306021
C5 0.0000000000 0.0000000000 -6.9242601761
C6 0.0000000000 2.4377336900 5.0640991723
End of basis
Basis set
H.ANO-S...2s.
H1 0.0000000000 4.0043085530 -3.8534714086
H2 0.0000000000 4.0326542950 0.7215314260
H3 0.0000000000 2.1467175630 7.0879851846
H4 1.5779129980 3.6622699270 4.5104123361
End of basis
RF-Input
reaction field
38.8 13.8 4
End of RF-Input
Well Int
4
1.0 5.0 15.8
1.0 3.5 16.8
1.0 2.0 18.8
1.0 1.4 20.8
End of Input
&SCF &END
TITLE
DMABN molecule
OCCUPIED
20 2 12 5
ITERATIONS
50
END OF INPUT
&RASSCF &END
TITLE
p-DMABN
SYMMETRY
1
SPIN
1
NACTEL
10 0 0
FROZEN
8 0 3 0
INACTIVE
12 1 9 1
RAS2
0 2 0 7
THRS
1.0E-06,1.0E-03,1.0E-03
ITER
50,25
LUMORB
HFRCtfld
1 1 1.796
END OF INPUT
&CASPT2 &END
Title
Solvent
Maxit
20
Lroot
1
RFPert
End of Input
&RASSCF &END
TITLE
p-DMABN
SYMMETRY
1
SPIN
1
NACTEL
10 0 0
FROZEN
8 0 3 0
INACTIVE
12 1 9 1
RAS2
0 2 0 7
THRS
1.0E-06,1.0E-03,1.0E-03
CIROot
1 2
2
Iterations
50,25
LumOrb
Hfrctfld
0 2 1.796
End of Input
&CASPT2 &END
Title
Solvent
Maxit
20
Lroot
2
RFPert
End of Input
This sequence of calculations is correct because we perform the CASPT2 calculations immediately after the RASSCF calculations for each root. They will not be correct if we first compute the two RASSCF states and then the two CASPT2 corrections, because in the COMFILE file the reaction field stored would be that from the last RASSCF calculation. The safer way is to save explicitly the COMFILE file each time.
One additional problem that can occur in excited state calculations
is that the RASSCF program does not converge easily for single root
calculations. One may have to do a state-average CASSCF calculation.
The program can compute the reaction field in average calculations
because at each iteration it takes the density matrix of the root
we specify under the keyword Hfrctfld to compute the reaction
field, althougth the optimized density matrix is the averaged one.
As this is only an approximation, the recommended procedure is to
increase as much as possible the weight of the computed state
in the average procedure. For instance, for the 2
A
state
in DMABN the RASSCF input could be:
&RASSCF &END
TITLE
p-DMABN molecule. 21A1 averaged state.
SYMMETRY
1
SPIN
1
NACTEL
10 0 0
FROZEN
8 0 3 0
INACTIVE
12 1 9 1
RAS2
0 2 0 7
CIROOT
2 2
1 2
1 9
LEVSHFT
0.5
ITER
50,25
CIMX
25
LUMORB
HFRctFld
0 2 1.801
END OF INPUT
The experimental and theoretical excitation energies, and computed
dipole moments for the ground state (GS) and 
2A
state of DMABN in different media are shown in Table 4.16
(see ref. [48]):
| Solvent | PT2 (eV) | Exp (eV) |  (D) |
 (D) |
|
| Gas phase | 4.41 | - | 6.6 | 14.2 | |
| Cyclohexane (=4.30) |
4.38 | 4.4 | 6.7 | 14.6 | |
| Butylchloride (=9.65) |
4.32 | 4.3 | 6.9 | 15.1 | |
| Acetonitrile (=38.8) |
4.31 | 4.2 | 7.0 | 15.2 | |
The large radius of the cavity and the empty space within the cavity lead to a clear underestimation of the solvent effects. In this case however the model accounts for the most important aspects of the interaction.
We now consider a difficult case. We want to compute both 
and
excited states in DMABN and we want to use different active
spaces in each case. For the states (2A and 1B
)
we use a (0207) space (
) with 10 active electrons and for the
states (1A and 1B)
a (1207) active space with 12 active electrons. This means we need
two different ground state calculations. In principle we could re-optimize
the cavity size for the new ground state but this is going to be a very
minor effect. The following is a Korn shell script designed to do the
explained calculations.
#!/bin/ksh ################################################################################ export Project=DMABN export Solvent=ACN export HomeDir=/u/$LOGNAME/$Project export TempDir=/temp/$LOGNAME export WorkDir=$TempDir/$RANDOM mkdir $WorkDir cd $WorkDir ln -fs $TempDir/$Project.OrdInt ORDINT ln -fs $TempDir/$Project.OneInt.$Solvent ONEINT #------------------------------------------------------------------------------# # Compute integrals # #------------------------------------------------------------------------------# molcas run seward $HomeDir/seward.$Project.$Solvent.input #------------------------------------------------------------------------------# # Compute RHF-SCF wavefunction # #------------------------------------------------------------------------------# ln -fs $TempDir/$Project.$Solvent.ScfOrb SCFORB molcas run scf $HomeDir/scf.$Project.input rm SCFORB #------------------------------------------------------------------------------# # Compute CASSCF state functions and CASPT2 energy corrections for the state # #------------------------------------------------------------------------------# # #------------------------------------------------------------------------------# # Ground State CASSCF and CASPT2 calculations # #------------------------------------------------------------------------------# Name_list='11A1_pipi 11A1_npi' for Name in $Name_list; do ln -fs $TempDir/$Project.$Solvent.ScfOrb INPORB ln -fs $TempDir/$Project.$Name.$Solvent.JobIph JOBIPH ln -fs $TempDir/$Project.$Name.$Solvent.RasOrb RASORB ln -fs $TempDir/$Project.$Name.$Solvent.ComFile COMFILE molcas run rasscf $HomeDir/rasscf.$Name.$Solvent.input molcas run rasread $HomeDir/rasread.input molcas run caspt2 $HomeDir/caspt2.rfpert.1.input rm RASORB rm INPORB rm JOBIPH rm COMFILE done #------------------------------------------------------------------------------# # Excited States CASSCF calculations # #------------------------------------------------------------------------------# Name_list='21A1 11B2' for Name in $Name_list; do #----------------------------------------------------------------# # Changing the name of the COMFILE # #----------------------------------------------------------------# cp $Project.11A1_pipi.$Solvent.ComFile $Project.$Name.$Solvent.ComFile #----------------------------------------------------------------# ln -fs $TempDir/$Project.11A1.$Solvent.RasOrb INPORB ln -fs $TempDir/$Project.$Name.$Solvent.JobIph JOBIPH ln -fs $TempDir/$Project.$Name.$Solvent.RasOrb RASORB ln -fs $TempDir/$Project.$Name.$Solvent.ComFile COMFILE molcas run rasscf $HomeDir/rasscfs.$Name.$Solvent.input molcas run rasread $HomeDir/rasread.input rm RASORB rm INPORB rm JOBIPH rm COMFILE done Name_list='11A2 11B1' for Name in $Name_list; do #----------------------------------------------------------------# # Changing the name of the COMFILE # #----------------------------------------------------------------# cp $Project.11A1_npi.$Solvent.ComFile $Project.$Name.$Solvent.ComFile #----------------------------------------------------------------# ln -fs $TempDir/$Project.11A1.$Solvent.RasOrb INPORB ln -fs $TempDir/$Project.$Name.$Solvent.JobIph JOBIPH ln -fs $TempDir/$Project.$Name.$Solvent.RasOrb RASORB ln -fs $TempDir/$Project.$Name.$Solvent.ComFile COMFILE molcas run rasscf $HomeDir/rasscfs.$Name.$Solvent.input molcas run rasread $HomeDir/rasread.input rm RASORB rm INPORB rm JOBIPH rm COMFILE done #------------------------------------------------------------------------------# # Excited States CASPT2 calculations # #------------------------------------------------------------------------------# Name_list='21A1' for Name in $Name_list; do ln -fs $TempDir/$Project.$Name.$Solvent.JobIph JOBIPH ln -fs $TempDir/$Project.$Name.$Solvent.ComFile COMFILE molcas run caspt2 $HomeDir/caspt2.rfpert.2.input rm JOBIPH rm COMFILE done Name_list='11B2' for Name in $Name_list; do ln -fs $TempDir/$Project.$Name.$Solvent.JobIph JOBIPH ln -fs $TempDir/$Project.$Name.$Solvent.ComFile COMFILE molcas run caspt2 $HomeDir/caspt2.rfpert.1.input rm JOBIPH rm COMFILE done Name_list='11A2 11B1' for Name in $Name_list; do ln -fs $TempDir/$Project.$Name.$Solvent.JobIph JOBIPH ln -fs $TempDir/$Project.$Name.$Solvent.ComFile COMFILE molcas run caspt2 $HomeDir/caspt2.rfpert.1.input rm JOBIPH rm COMFILE done #------------------------------------------------------------------------------# cd $TempDir rm -r $WorkDir exit
Notice that the two-electron integral file is the same independent of the cavity size, well integrals used or translational movement of the molecule. Therefore, if the well parameters are changed only the one-electron integral file ONEINT need to be recomputed using the option ONEOnly in the SEWARD input. In the previous script we have named the ONEINT as $Project.OneInt.$Solvent, but not because it depends on the solvent but on the cavity radius which should be different for each solvent when the well potential is used.