$TRANST group (relevant for RUNTYP=TRANSITN)
(only for CITYP=GUGA or MPLEVL=2)
This group controls the evaluation of the radiative
transition moment, or spin orbit coupling (SOC). An SOC
calculation can be based on variational CI wavefunctions,
using GUGA CSFs, or based on 2nd order perturbation theory
using the MCQDPT multireference perturbation theory.
These are termed SO-CI and SO-MCQDPT below. The orbitals
are typically obtained by MCSCF computations, and since
the CI or MCQDPT wavefunctions are based on those MCSCF
states, the zero-th order states are referred to below as
the CAS-CI states. SOC jobs prepare a model Hamiltonian
in the CAS-CI basis, and diagonalize it to produce spin-
mixed states, which are linear combinations of the CAS-CI
states. If scalar relativistic corrections were included
in the underlying spin-free wavefunctions, it is possible
either to include or to neglect similar corrections to the
spin-orbit integrals, see keyword NESOC in $RELWFN.
An input file to perform SO-CI will contain
SCFTYP=NONE CITYP=GUGA MPLEVL=0 RUNTYP=TRANSITN
while a SO-MCQDPT calculation will have
SCFTYP=NONE CITYP=NONE MPLEVL=2 RUNTYP=TRANSITN
The SOC job will compute a Hamiltonian matrix as the sum
of spin-free terms and spin-orbit terms, H = H-sf + H-so.
For SO-CI, the matrix H-sf is diagonal in the CAS-CI state
basis, with the LS-coupled CAS-CI energies as the diagonal
elements, and H-so contains only off-diagonal couplings
between these LS states,
H-sf = CAS-CI spin-free E
H-so = CAS SOC Hamiltonian (e.g. HSO1, HSO2P, HSO2)
For SO-MCQDPT, the additional input PARMP defines these
matrices differently. For PARMP=0, the spin-free term
has diagonal and off-diagonal MCQDPT perturbations:
H-sf - CAS-CI spin-free E + 2nd order spin-free MCQDPT
H-so - CAS SOC Hamiltonian
For PARMP not equal to 0, the spin orbit operator is also
included into the perturbing Hamiltonian of the MCQDPT:
H-sf - CAS-CI spin-free E + 2nd order spin-free MCQDPT
H-so - CAS SOC Hamiltonian + 2nd order SO-MCQDPT
Pure transition moment calculations (OPERAT=DM) are
presently limited to CI wavefunctions, so please use only
CITYP=GUGA MPLEVL=0. The transition moments computed by
SO-MCQDPT runs (see TMOMNT flag) will form the transition
density for the CAS-CI zeroth order states rather than the
1st order perturbed wavefunctions.
Please see REFS.DOC for additional information on what
is actually a fairly complex input file to prepare.
OPERAT selects the type of transition being computed.
= DM calculates radiative transition moment
between states of same spin, using
the dipole moment operator. (default)
= HSO1 one-electron Spin-Orbit Coupling (SOC)
= HSO2P partial two electron and full 1e- SOC,
namely core-active 2e- contributions are
computed, but active-active 2e- terms
are ignored. This generally captures
>90% of the full HSO2 computation, but
with spin-orbit matrix element time
similar to the HSO1 calculation.
= HSO2 one and two-electron SOC, this is the
full Pauli-Breit operator.
= HSO2FF one and two-electron SOC, the form factor
method gives the same result as HSO2, but
is more efficient in the case of small
active spaces, small numbers of CAS-CI
states, and large atomic basis sets.
This final option applies only to SO-CI.
PARMP = controls inclusion of the SOC terms in SO-MCQDPT,
for OPERAT=HSO1 (default=1) or for HSO2P/HSO2
(default=3) only.
0 - no SOC terms should be included in the MCQDPT
corrections at 2nd order, but they will be
included in the CAS states on which the MCQDPT
(i.e. up to 1st order)
1 - include the 1e- SOC perturbation in MCQDPT
-1 - defined under "3", read on...
3 - full 1-electron and partial 2-electron in the
form of the mean field perturbation (this is
very similar to HSO2P, but in the MCQDPT2
perturbation). Only doubly occupied orbitals
(NMODOC) are used for the core 2e terms.
If the option is set to -1, then all core
orbitals (NMOFZC+NMODOC) are used. Neither
calculation includes extra diagrams including
filled orbitals, so both are "partial".
PARMP=3 (or -1) has almost no extra cost compared to
PARMP=1, but can only be used with OPERAT=HSO2 or HSO2P.
The options -1 and 3 are not rigorously justified, contrary
to HOS2P for a SO-CI, as 2e integrals with 2 core indices
appear in the second order in two ways. There is a mean-
field addition to 1e integrals, which is included when you
choose PARMP=3 or -1. But, there are separate terms from
additional diagrams that are not implemented, so that there
is some imbalance in including the partial 2e correction.
Nevetheless, it may be better to include such "partial"
partial 2e contributions than not to. Note that at first
order in the energy (the CAS-CI states) the N-electron
terms are treated exactly as specified by OPERAT.
NFFBUF = sets buffer size for form factors in SO-MCQDPT.
(applies only to OPERAT = HSO1, HSO2 or HSO2P).
This is a very powerful option that speeds up
SO-MCQDPT calculations by precomputing the total
multiplicative factor in front of each diagram so
that the latter is computed only once (this is in
fact what happens in MCQDPT). It is not uncommon
for this option to speed up calculations by a
factor of 10. Since this option forces running
the SO-CASCI part twice (due to the SO-MCQDPT
Hamiltonian being non-Hermitian), it is possible
that in rare cases NFFBUF=0 may perform similarly
or better. The upper bound for NFFBUF is NACT**2,
where NACT=NOCC-NFZC. Due to the sparseness of
the coupling constants it is usually sufficient to
set NFFBUF to 3*NACT. To use the older way of
dynamically computing form factors and diagrams on
the fly, set NFFBUF to 0. Default: 3*(NOCC-NFZC)
It is advisable to tighten up the convergence criteria in
the $MCQDx groups since SOC is a fairly small effect, and
the spin-free energies should be accurately computed, for
example THRCON=1e-8 THRGEN=1e-10.
PARMP has a rather different meaning for OPERAT=HSO2FF:
It refers to the difference between ket and bra's Ms,
-1 do matrix elements for ms=-1 only
0 do matrix elements for ms=0 only
1 do matrix elements for ms=1 only
-2 do matrix elements for all ms (0, 1, and -1),
which is the default.
-3 calculates form factors so they can be saved
* * * next defines the orbitals and wavefunctions * * *
NUMCI = For SO-CI, this parameter tells how many CI
calculations to do, and therefore defines how
many $DRTx groups will be read in.
For SO-MCQDPT, this parameter tells how many
MCQDPT calculations to do, and therefore defines
how many $MCQDx groups will be read in.
(default=1)
IROOTS, IVEX, NSTATE, and ENGYST below will all
have NUMCI values. NUMCI may not exceed 64.
You may wish to define one $DRTx or $MCQDx group for each
spatial symmetry representation occurring within each spin
multiplicity, as the use of symmetry during these separate
calculations may make the entire job run much faster.
NUMVEC = the meaning is different depending on the run:
a) spin-orbit CI (SO-CI),
Gives the number of different MO sets. This can
be either 1 or 2, but 2 can be chosen only for
FORS/CASSCF or FCI wavefunctions. (default=1)
If you set NUMVEC=2 and you use symmetry in any
of the $DRTx groups, you may have to use STSYM
in the $DRTx groups since the order of orbitals
from the corresponding orbital transformation
is unpredictable.
b) spin-orbit perturbation (SO-MCQDPT),
The option to have different MOs for different
states is not implemented, so your job will have
only one $VEC1 group, and IVEX will not normally
be input. The absolute value of NUMVEC should be
be equal to the value of NUMCI above. If NUMVEC
positive, the orbitals in the $VEC1 will be used
exactly as given, whereas if NUMVEC is a negative
number, the orbitals will be canonicalized
according to IFORB in $MCQDx. Using NUMVEC=-NUMCI
and IFORB=3 in all $MCQDx to canonicalize over all
states is recommended.
Note that $GUESS is not read by this RUNTYP! Orbitals must
be in $VEC1 and possibly $VEC2 input groups.
NFZC = For SO-CI, this is equal to NFZC in each $DRTx
group. When NUMVEC=2, this is also the number of
identical core orbitals in the two vector sets.
For SO-MCQDPT, this should be NMOFZC+NMODOC given
in each of the $MCQDx groups.
The default is the number of AOs given in $DATA,
this is not very reasonable.
NOCC = the number of occupied orbitals. For SO-CI this
should be NFZC+NDOC+NALP+NAOS+NBOS+NVAL, but
add the external orbitals if the CAS-CI states
are CI-SD or FOCI or SOCI type instead of CAS.
For SO-MCQDPT enter NUMFZC+NUMDOC+NUMACT.
The default is the number of AOs given in $DATA,
which is not usually correct.
Note: IROOTS, NSTATE, ENGYST, IVEX contain NUMCI values.
IROOTS = array containing the number of CAS-CI states to
be used from each CI or MCQDPT calculation.
The default is 1 for every calculation, which is
probably not a correct choice for OPERAT=DM runs,
but is quite reasonable for the HSO operators.
The total number of states included in the SOC
Hamiltonian is the summation of the NUMCI values
of IROOTS times the multiplicity of each CI or
MCQDPT. See also ETOL/UPPREN.
NSTATE = array containing the number of CAS-CI states to be
found by diagonalising the spin-free Hamiltonians.
Of these, the first IROOTS(i) states will be used
to find transition moments or SOC. Obviously,
enter NSTATE(i) >= IROOTS(i).
The default for NSTATE(i) is IROOTS(i), but might
be bigger if you are curious about the additional
energies, or to help the Davidson diagonalizer.
NSTATE is ignored by SO-MCQDPT runs, and you must
ensure that your IROOTS input corresponds to the
KSTATE option in $MCQDx.
ETOL = energy tolerance for CI state elimination.
This applies only to SO-CI and OPERAT=HSO1,2,2P.
After each CI finds NSTATE(i) CI roots for each
$DRTx, the number of states kept in the run is
normally IROOTS(i), but ETOL applies the further
constraint that the states kept be within ETOL of
the lowest energy found for any of the $DRTx.
The default is 100.0 Hartree, so that IROOTS is
the only limitation.
UPPREN = similar to ETOL, except it is an absolute energy,
instead of an energy difference.
IVEX = Array of indices of $VECx groups to be used for
each CI calculation. The default for NUMVEC=2 is
IVEX(1)=1,2,1,1,1,1,1..., and of course for
NUMVEC=1, it is IVEX(1)=1,1,1,1,1...
This applies only to CITYP=GUGA jobs.
ENGYST = energy values to replace the spin-free energies.
This parameter applies to SO-CI only.
A possible use for this is to use first or second
order CI energies (FOCI or SOCI in $DRT) on the
diagonal of the Hamiltonian (obtained in some
earlier runs) but to use only CAS wavefunctions
to evaluate off diagonal HSO matrix elements.
The CAS-CI is still conducted to get CI coefs,
needed to evaluate the off diagonal elements.
Enter MXRT*NUMCI values as a square array, by the
usual FORTRAN convention (that is, MXRT roots of
$DRT1, MXRT roots of $DRT2 etc), in hartrees, with
zeros added to fill each column to MXRT values.
MXRT is the maximum value in the IROOTS array.
(the default is the computed CAS-CI energies)
See B.Schimmelpfennig, L.Maron, U.Wahlgren,
C.Teichteil, H.Fagerli, O.Gropen Chem.Phys.Lett.
286, 261-266(1998).
* * * the next pertain only to spin-orbit runs * * *
ISTNO if given as positive values:
an array of one or two state indices which govern
computation of the density matrix of one state,
or the transition density of two states.
if given as negative values:
one state-averaged density with equal weights.
ISTNO(1)=5 state-specific density of state 5
ISTNO(1)=1,2 transition density between 1 and 2
ISTNO(1)=-1,-6 state-average all states 1 to 6
The default is ISTNO(1)=0,0 meaning no density.
Computation of the density gives access to the
full Gaussian property package, except Mulliken
populations. At present, computation of the
transition density does just that, without any
oscillator strengths. If the computation is of
SO-MCQDPT type, the density or transition density
that is computed will be that for the unperturbed
SO-CASCI states.
DEGTOL = array of two tolerances to help define what states
are considered degenerate. This is ignored except
for linear molecules or atoms. The purpose is to
decide what states are grouped together during the
determination of simultaneous eigenstates of the
spin-orbit Hamiltonian and Jz. DEGTOL(1) is in
wavenumbers, and defines which spin-orbit states
have the same energy. DEGTOL(2) is in units of
electrons, and defines which natural orbitals are
considered to be degenerate. If the Jz values in
your run seem incorrect, tighten or relax the two
degeneracy tolerances to get the correct groupings
of the states. Default= 0.02,0.002
RSTATE = sets the zero energy level
format: ndrt*1000+iroot for adiabatic state (root)
0000 sets zero energy to the lowest diabatic root
default: 1001 (1st root in $DRT1 or $MCQD1)
ZEFTYP specifies effective nuclear charges to use.
= TRUE uses true nuclear charge of each atom,
except protons are removed if an ECP basis
is being used (default).
= 3-21G selects values optimized for the 3-21G
basis, but these are probably appropriate
for any all electron basis set. Rare gases,
transition metals, and Z>54 will use the
true nuclear charges.
= SBKJC selects a set obtained for the SBKJC ECP
basis set, specifically. It may not be
sensible to use these for other ECP sets.
Rare gases, lanthanides, and Z>86 will use
the true nuclear charges.
ZEFF = an array of effective nuclear charges, overriding
the charges chosen in ZEFTYP.
Note that effective nuclear charges can be used for
any HSO type OPERAT, but traditionally these are used
mainly for HSO1 as an empirical correction to the
omission of the 2e- term, or to compensate for missing
core orbitals in ECP runs.
ONECNT = uses a one-center approximation for SOC integrals:
= 0 compute all SOC integrals without approximations
= 1 compute only one-center 1e and 2e SOC integrals
= 2 compute all 1e, but only one-center 2e integrals
Numerical tests indicate the error of the one-center
approximation (ONECNT=1) is usually on the order of a
few wavenumbers for Li-Ne (a bit larger for F) and its
errors appear to become negligible for anything heavier
than Ne. ONECNT=1 appears to give a better balanced
description than ONECNT=2. Very careful users can check
how well the approximation works for their particular
system by using ONECNT=0, then ONECNT=1, to compare
the results. One important advantage of ONECNT=1/2 is
that this removes the dependence of SOC 2e integrals
upon the molecular geometry. This means the program
needs to compute SOC 2e integrals only once for a given
set of atoms and then they can be read by using SOC
integral restart. RUNTYP=SURFACE automatically takes
advantage of this fact.
JZ controls the calculation of Jz eigenvalues
= 0 do not perform the calculation
= 1 do the calculation
By default, Jz is set to 1 for molecules that are
recognised as linear (this includes atoms!).
Jz cannot be computed for nonlinear molecules.
The matrix of Jz=Lz+Sz operator is constructed
between spin-mixed states (eigenvalues of Hso).
Setting Jz to 1 can enforce otherwise avoided (by
symmetry) calculations of SOC matrix elements.
JZ applies only to HSO1,2,2P.
TMOMNT = flag to control computation of the transition
dipole moment between spin-mixed wavefunctions
(that is, between eigenvectors of the Pauli-Breit
Hamiltonian). Applies only to HSO1,2,2P.
(default is .FALSE.)
SKIPDM = flag to omit(.TRUE.) or include(.FALSE.) dipole
moment matrix elements during spin-orbit coupling.
Usually it takes almost no addition effort to
calculate excluding some cases when the
calculation of forbidden by symmetry spin-orbit
coupling matrix elements may have to be
performed since and are computed
simultaneously. Applies only to HSO1,2,2P.
Since the lack of a MCQDPT density matrix means
there are no MCQDPT dipole moments at present,
SO-MCQDPT jobs will compute the dipole matrix
elements for the CAS-CI states only. However,
the dipole moments in the spin-mixed states will
be computed with the MCQDPT mixing coefficients.
(default is .TRUE.)
IPRHSO = controls output style for matrix elements (HSO*)
=-1 do not output individual matrix elements
otherwise these are accumulative:
= 0 term-symbol like kind of labelling:
labels contain full symmetry info (default)
= 1 all states are numbered consequently within each
spin multiplicity (ye olde style)
= 2 output only nonzero (>=1e-4) matrix elements
PRTPRM = flag to provide detailed information about the
composition of the spin-mixed states in terms of
adiabatic states. This flag also provides similar
information about Jz (if JZ set).
(default is .FALSE.)
LVAL = additional angular momentum symmetry values:
For the case of running an atom:
LVAL is an array of the L values (L**2 = L(L+1))
for each $MCQDx/$DRTx (L=0 is S, 1 is P, etc.)
For the case of running a linear molecule:
LVAL is an array that gives the |Lz| values. Note
that real-valued wavefunctions (e.g. Pi-x, Pi-y)
have Lz and -Lz components mixed, so you should
input |Lz| as 1 and 1 for both Pi-x and Pi-y.
This parameter should not be given for a nonlinear
polyatomic system.
Default: all set to -1 (that is, do not use these
additional symmetry labels. It is the user's
responsibility to ensure the values' correctness.
Note that for SO-MCQDPT useful options in $MCQDPT
are NDIAOP and KSTATE. They enable efficient
separation of atomic/linear symmetry irreps).
It is acceptable to set only some values and leave
others as -1, if only some groups have definite
values. Note that normally Lz values are printed
at the end of the log file, so its easy to double
check the initial values for LVAL. For the case
of atoms LVAL drastically reduces the CPU time, as
it reduces a square matrix to tridiagonal form.
For the case of linear molecules the savings at
the spin-orbit level are somewhat less, but they
are usually quite significant at the preceding
spin-free MCQDPT step.
MCP2E = Model Core Potential SOC 2e contributions.
Note that MCP 1e contributions are handled as in
case of all-electron runs because MCP orbitals
contain all proper nodes).
= 0 do not add the MCP 2e core-active contribution,
but add any other 2e- terms asked for by OPERAT.
= 1 add this contribution, but no other 2e SOC term.
This is recommended, and the default.
= 2 add this contribution and the 2e- contributions
requested by OPERAT, for any e- which are being
treated by quantum mechanics (not MCP cores).
Note that for MCP2E=0 and 2, HSO2, HSO1, HSO2P
values of OPERAT are supported for the explicit
2e- contributions. The recommended approach is to
assume that MCP alone can capture all the 2e SOC,
for this use MCP2E=1 OPERAT=HSO2P. The entire 2e-
contribution is achieved with MCP2E=2 OPERAT=HSO2.
If your MCP leaves out many core electrons as
particles, MCP2E=2 OPERAT=HSO2P can be tested to
see if it adds a sizable amount to SOC, compared
to MCP2E=1 OPERAT=HSO2P).
MCP2E=2 OPERAT=HSO1 is an illegal combination.
MCP2E=1 OPERAT=HSO1 is illogical since the MCP 2e
integrals are computed but not used anywhere.
The following table explains MCP2E and gives all
useful combinations:
MCP2E/OPERAT 2e SOC contributions SOC 2e ints
2 HSO2 MCPcore-CIact + CIcore-CIact MCP+basis
+ CIact-CIact
2 HSO2P MCPcore-CIact + CIcore-CIact MCP+basis
1 HSO2P MCPcore-CIact MCP
using the following orbital space definitions:
MCPcore orbitals whose e- are replaced by MCP
CIcore always doubly occupied
CIact MOs allowed to have variable occupation
* * * expert mode HSO control options * * *
MODPAR = parallel options, which are independent bit
options, 0=off, 1=on. Bit 1 refers only to
HSO2FF, bit 2 to HSO1,2,2P. Enter a decimal
value 0, 1, 2, 3 meaning binary 00, 01, 10, 11.
bit 1 = 0/1 (HSO2FF) uses static/dynamic load balancing in
parallel if available, otherwise use static
load balancing. Dynamic algorithm is usually
faster but may utilize memory less efficiently,
and I/O can slow it down. Also, dynamical
algorithm forces SAVDSK=.F. since its
unique distribution of FFs among nodes implies
no savings from precalculating form factors.
bit 2 = 0/1 (HSO1,2,2P) duplicate/distribute SOC integrals
in parallel. If set, 2e AO integrals and the
four-index transformation are divided over
nodes (distributed), and SOC MO integrals are
then summed over nodes.
The default is 3, meaning both bits are set on (11)
PHYSRC = flag to force the size of the physical record to
be equal to the size of the sorting buffers.
This option can have a dramatic effect on the
efficiency. Usually, setting PHYSRC=.TRUE. helps
if the code complains that low memory enforces
SLOWFF=.TRUE., or you set it yourself. For large
active spaces and large memory (more precisely, if
RECLEN is larger than the physical record size)
PHYSRC=.TRUE. can slow the code down. Setting
PHYSRC to .true. forces SLOWFF to be .false.
See MODPAR. (default .FALSE.) (only with HSO2FF)
RECLEN = specifies the size of the record on file 40,
where form factors are stored. This parameter
significantly affects performance.
If not specified, RECLEN have to be guessed,
and the guess will usually be either an
overestimate or underestimate. If the former
you waste disk space, if the latter the program
aborts. Note that RECLEN will be different for
each pair of multiplicities and you must specify
the maximum for all pairs. The meaning of this
number is how many non-zero form factors are
present given four MO indices. You can decrease
RECLEN if you are getting a message "predicted
sorting buffer length is greater than needed..."
Default depends on active space. (only HSO2FF)
SAVDSK = flag to repeat the form factor calculation twice.
This avoids wasting disk space as the actually
required record size is found during the 1st run.
(default=.FALSE.) (only with HSO2FF)
SLOWFF = flag to choose a slower FF write-out method.
By default .FALSE., but this is turned on if:
1) not enough memory for the fast way is available
2) the maximum usable memory is available, as when
the buffer is as large as the maximum needed,
then the "slow FF" algorythm is faster.
Generally SLOWFF=.true. saves up to 50% or so of
disk space. See PHYSRC. (only with HSO2FF)
ACTION controls disk file DAFL30 reuse.
= NORMAL calculate the form factors in this run.
= SAVE calculate, and store the form factors on
disk for future runs with the same active
space characteristics.
= READ read the form factors from disk from an
earlier run which used SAVE.
(default=NORMAL) (only with HSO2FF)
Note that currently in order to use ACTION =
SAVE or READ you should specify MS= -1, 0, or 1
* * * some control tolerances * * *
NOSYM= -1 forces use of symmetry-contaminated orbitals
symmetry analysis, otherwise the same as NOSYM=0
= 0 fully use symmetry
= 1 do not use point group symmetry, but still use
other symmetries (Hermiticity, spin).
= 2 use no symmetry. Also, include all CSFs for
HSO1, 2, 2P.
= 3 force the code to assume the symmetry specified
in $DATA is the same as in all $DRTx groups, but
is otherwise identical to NOSYM=-1. This option
saves CPU time and money(memory). Since the $DRT
works by mapping non-Abelian groups into their
highest Abelian subgroup, do not use NOSYM=3 for
non-Abelian groups.
SYMTOL = relative error for the matrix elements. This
parameter has a great impact upon CPU time, and
the default has been chosen to obtain nearly
full accuracy while still getting good speedups.
(default=1.0E-4)
* * * the remaining parameters are not so important * * *
PRTCMO = flag to control printout of the corresponding
orbitals. (default is .FALSE.)
HSOTOL = HSO matrix elements are considered zero if they
are smaller than HSOTOL. This parameter is used
only for print-out and statistics.
(default=1.0E-1 cm-1)
TOLZ = MO coefficient zero tolerance (as for $GUESS).
(default=1.0E-8)
TOLE = MO coefficient equating tolerance (as for
$GUESS). (default=1.0E-5)
==========================================================
* * * * * * * * * * * * * * * * * * *
For information on RUNTYP=TRANSITN,
see the 'further information' section
* * * * * * * * * * * * * * * * * * *
587 lines are written.
Edited by Shiro KOSEKI on Thu Mar 5 10:25:38 2020.