SEWARD usually runs after program GATEWAY. In the same time, any input used
in GATEWAY can be placed into SEWARD input. However, it is recommended to
specify all information about the molecule and the basis set in GATEWAY input.

Label  Explanation 
'Mltpl nn'  the nn'th order Cartesian multipole moments. 
'MltplS '  the overlap matrix used in the semiempirical NDDO method. 
'Kinetic '  the kinetic energy integrals. 
'Attract '  the electron attraction integrals. 
'AttractS'  the electron attraction integrals used in the semiempirical NDDO method. 
'PrjInt '  the projection integrals used in ECP calculations. 
'M1Int '  the M1 integrals used in ECP calculations. 
'M2Int '  the M2 integrals used in ECP calculations. 
'SROInt '  the spectrally resolved operator integrals used in ECP calculations. 
'XFdInt '  the external electric field integrals. 
'MassVel '  the massvelocity integrals. 
'Darwin '  the Darwin oneelectron contact integrals. 
'Velocity'  the velocity integrals. 
'EF0nnnnn'  the electric potential at center nnnnn. 
'EF1nnnnn'  the electric field at center nnnnn. 
'EF2nnnnn'  the electric field gradient at center nnnnn. 
'AngMom '  the angular momentum integrals. 
'DMS '  the diamagnetic shielding integrals. 
'Wellnnnn'  the nnnn'th set of spherical well integrals. 
'OneHam '  the oneelectron Hamiltonian. 
'AMProd '  the hermitized product of angular momentum integrals. 
'AMFI '  the atomic mean field integrals. 
The input for each module is preceded by its name like:
&SEWARD
Argument(s) to a keyword, either individual or composed by several entries, can be placed in a separated line or in the same line separated by a semicolon. If in the same line, the first argument requires an equal sign after the name of the keyword.
Keyword  Meaning 
TITLe  One line of title card follows. 
TEST  SEWARD will only process the input and generate a nonzero return code. 
ONEOnly  SEWARD will not compute the twoelectron integrals. 
NODKroll  SEWARD will not compute DouglasKroll integrals. 
DIREct  Prepares for later integraldirect calculations. As with keyword OneOnly, SEWARD will evaluate no twoelectron integrals. 
EXPErt  Sets ``expert mode'', in which various default settings are altered. Integraldirect calculations will be carried out if the twoelectron integral file is unavailable. 
CHOLesky  SEWARD will Cholesky decompose the twoelectron integrals using default configuration (in particular, the decomposition threshold is 1.0d4) of the decomposition driver. The configuration may be tailored using the ChoInput section. Default is to not decompose. 
LOW Cholesky  SEWARD will Cholesky decompose the twoelectron integrals using lowaccuracy (threshold 1.0d4) configuration of the decomposition driver (the configuration may be tailored using the ChoInput section). Default is to not decompose. 
MEDIum Cholesky  SEWARD will Cholesky decompose the twoelectron integrals using mediumaccuracy (threshold 1.0d6) configuration of the decomposition driver (the configuration may be tailored using the ChoInput section). Same as Cholesky. Default is to not decompose. 
HIGH Cholesky  SEWARD will Cholesky decompose the twoelectron integrals using highaccuracy (threshold 1.0d8) configuration of the decomposition driver (the configuration may be tailored using the ChoInput section). Default is to not decompose. 
FAKE CD/RI  If CD/RI vectors are already available, SEWARD will not redo work! 
JMAX  The integer entry on the next line is the highest rotational quantum number for which SEWARD will compute the rotational energy within the rigid rotor model. The default value is 5. 
SYMMetry  See the the description in the manual for the program GATEWAY 
BASIs Set  See the the description in the manual for the program GATEWAY 
ZMAT  See the the description in the manual for the program GATEWAY 
NOGUessorb  Disable automatic generation of starting orbitals with the GuessOrb procedure. 
NODElete  Do not delete any orbitals automatically. 
SDELete  Set the threshold for deleting orbitals based on the eigenvalues of the overlap matrix. All eigenvalues with eigenvectors below this threshold will be deleted. If you want no orbitals deleted use keyword NODElete. 
TDELete  Set the threshold for deleting orbitals based on the eigenvalues of the kinetic energy matrix. All eigenvalues with eigenvectors above this threshold will be deleted. If you want no orbitals deleted use keyword NODElete. 
ECPShow  Force SEWARD to print ECP parameters. 
AUXShow  Force SEWARD to print auxiliary basis set parameters. 
BSSHow  Force SEWARD to print basis set parameters. 
VERbose  Force SEWARD to print a bit more verbose. 
Keyword  Meaning 
MULTipoles  Entry which specifies the highest order of the multipole for which integrals will be generated. The default center for the dipole moment operator is the origin. The default center for the higher order operators is the center of the nuclear mass. The default is to do up to quadrupole moment integrals (2). 
CENTer  This option is used to override the default selection of the origin of the multipole moment operators. On the first entry add an integer entry specifying the number of multipole moment operators for which the origin of expansion will be defined. Following this, one entry for each operator, the order of the multipole operator and the coordinates of the center (in au) of expansion are specified. 
SDIPole  Supplement ONEINT for transition dipole moment calculations, i.e. dipole moment and velocity integrals will be computed. This option should be used whenever the RASSI program is used to compute transition moments, so that the transition moments can be evaluated in both velocity and length representation. 
ANGM  Supplement ONEINT for transition angular momentum calculations. Entry which specifies the angular momentum origin (in au). 
DSHD  Requests the computation of diamagnetic shielding integrals. The first entry specifies the gauge origin. Then follows an integer specifying the number of points at which the diamagnetic shielding will be computed. If this entry is zero, the diamagnetic shielding will be computed at each nucleus. If nonzero, then the coordinates (in au) for each origin has to be supplied, one entry for each origin. 
RELInt  Requests the computation of massvelocity and oneelectron Darwin contact term integrals for the calculation of a first order correction of the energy with respect to relativistic effects. 
AMPR  Request the computation of angular momentum product integrals. The keyword is followed by values which specifies the angular momentum origin (in au). 
RXXPyy  Request arbitrary scalar relativistic DouglasKrollHess (DKH) correction to the oneelectron Hamiltonian
and the socalled picturechange correction to the property integrals (multipole moments
and electronic potential related properties).
Here XX represents the order of the DKH correction to the oneelectron Hamiltonian and
yy the order of the picturechange correction. The character P denotes the parameterization
used in the DKH procedure.
The possible parametrizations P of the unitary transformation used
in the DKH transformation supported by MOLCAS are:

RX2C  Request the scalar relativistic X2C (eXacttwoComponent) corrections to the oneelectron Hamiltonian as well as the property integrals. 
RBSS  Request the scalar relativistic BSS (BaryszSadlejSnijders) corrections to the oneelectron Hamiltonian as well as the property integrals. The noniterative scheme is employed for the construction of BSS transformation. 
NOAMfi  Explicit request for no computation of atomic meanfield integrals. 
AMFI  Explicit request for the computation of atomic meanfield integrals (used in subsequent spinorbit calculations). These integrals are computed by default for the ANORCC and ANODK3 basis sets. 
EPOT  An integer follows which represents the number of points for which the electric potential will be computed. If this number is zero, the electric field acting on each nucleus will be computed. If nonzero, then the coordinates (in au) for each point have to be supplied, one entry for each point. 
EFLD  An integer follows which represents the number of points for which the electric potential and electric field will be computed. If this number is zero, the electric field acting on each nucleus will be computed. If nonzero, then the coordinates (in au) for each point have to be supplied, one entry for each point. 
FLDG  An integer required which represents the number of points for which the electric potential, electric field and electric field gradient will be computed. If this number is zero, the electric field gradient acting on each nucleus will be computed. If nonzero, then the either coordinates (in au) for each point or labels for each atom center have to be supplied, one entry for each point. In case a label is supplied it must match one of those given previous in the input during specification of the coordinates of the atom centers. Using a label instead of a coordinate can e.g. be useful in something like a geometry optimization where the coordinate isn't known when the input is written. 
Grid Input  Specification of numerical quadrature parameters, consult the numerical quadrature section of this manual. 
Keyword  Meaning 
VECTors  Requests a property calculation. For this purpose a file, INPORB, must be available, which contains the MO's and occupation numbers of a wave function. 
ORBCon  The keyword will force SEWARD to produce a list of the orbital contributions to the properties being computed. The default is to generate a compact list. 
THRS  The real entry on the following line specifies the threshold for the occupation number of an orbital in order for the OrbCon option to list the contribution of that orbital to a property. The default is 1.0d6. 
Keyword  Meaning 
NOPAck  The twoelectron integrals will not be packed. The default is to pack the twoelectron integrals. 
PKTHre  An entry specifies the desired accuracy for the packing algorithm, the default is 1.0d10. 
STDOut  Generate a twoelectron integral file according to the standard of version 1 of MOLCAS. The default is to generate the twoelectron integrals according to the standard used since version 2 of MOLCAS. 
THREshold  Threshold for writing integrals to disk follows. The default is 1.0d10. 
CUTOff  Threshold for ignoring the calculation of integrals based on the pair prefactor follows. The default is 1.0d10. 
This section contains keyword which control the radial numerical integration of the diffuse basis functions describing the scattered electrons in the variational Rmatrix approach. The activation of this option is controlled by that the center of the diffuse basis is assigned the unique atom label DBAS.
Keyword  Meaning 
RMAT  Radius of the Rmatrix sphere (in Bohr). This sphere is centered at the coordinate origin. The default is 10 Bohr. 
RMEA  Absolute precision in radial integration. The default is 1d9. 
RMER  Relative precision in radial integration. The default is 1d14. 
RMQC  Effective charge of the target molecule. This is the effective charge seen by the incident electron outside of the Rmatrix sphere. The default is 0d0. 
RMDI  Effective dipole of the target molecule. This is the effective dipole seen by the incident electron outside of the Rmatrix sphere. The default is (0d0,0d0,0d0). 
RMEQ  Minimal value of the effective charge of the target molecule to be considered. This is also the minimal value of the components of the effective dipole to be considered. Default is 1d8 
RMBP  Parameter used for test purposes in the definition of the Bloch term. Default is 0d0. 
CELL  Defines the three vectors of the unit cell (,,). The optional keyword Angstrom before the definition of vectors would read data in Å. Must consist of three entries (four in the case of Å) which correspond to coordinates of the vectors. All the atoms which are defined after that key are considered as the atoms of the cell. 
SPREad  Three integer numbers n_{1}, n_{2}, n_{3} which define the spread of the unit cell along the unit cell vectors. For example, would add all cell's atoms translated on , , , . This key must be placed before the definition of the unit cell atoms. 
Below follows an input for the calculation of integrals of a carbon atom. The comments in the input gives a brief explanation of the subsequent keywords.
&SEWARD
Title= This is a test deck!
* Remove integrals from a specific irreps
Skip= 0 0 0 0 1 1 1 1
* Requesting only overlap integrals.
Multipole= 0
* Request integrals for diamagnetic shielding
DSHD= 0.0 0.0 0.0; 1; 0.0 0.0 0.0
* Specify a title card
* Request only one}electron integrals to be computed
OneOnly
* Specify group generators
Symmetry= X Y Z
* Specify basis sets
Basis set
C.ANOL...6s5p3d2f.
Contaminant d
C 0.0 0.0 0.0
End of basis
The label, which defines the basis set for a given atom or set of atoms, is given as input after the keyword Basis set. It has the following general structure (notice that the last character is a period):
where the different identifiers have the following meaning:
Identifier  Meaning 
atom  Specification of the atom by its chemical symbol. 
type  Gives the type of basis set (ANO, STO, ECP, etc.) according to specifications given in the basis set library, vide supra. Observe that the upper cased character of the type label defines the file name in the basis directory. 
author  First author in the publication where that basis set appeared. 
primitive  Specification of the primitive set (e.g. 14s9p4d3f). 
contracted  Specification of the contracted set to be selected. Some basis sets allow only one type of contraction, others all types up to a maximum. The first basis functions for each angular momentum is then used. Note, for the basis set types ANO and ECP, onthefly decontraction of the most diffuse functions are performed in case the number of contracted functions specified in this field exceeds what formally is specified in the library file. 
aux  Specification of the type of AIMP, for instance, to choose between nonrelativistic and relativistic core AIMP's. 
atom
, type
,
and contracted
have to be included in
the label. The others can be left out. However, the periods have to be
kept. Example  the basis set label
`C.anos...4s3p2d.
'
will by MOLCAS be
interpreted as
`C.anos.Pierloot.10s6p3d.4s3p2d.
',
which is the first
basis set in the ANOS file in the
basis directory that fulfills the specifications given.
where Charge is the nuclear charge, lAng is the highest angular momentum quantum number, nPrim is the number of primitive functions (exponents) for a given shell, and nContr is the number of contracted functions for a given shell.
The following is an example of a DPZ basis set for carbon. Normally, however, the basis set will be read from a library file following the specified label (like, e.g., C.DZP...4s2p1d.), and not be inserted inline at the input file.
Basis set  Start defining a basis set
C.DZP.Someone.9s5p1d.4s2p1d. / inline  Definition in input stream
6.0 2  charge, max lquantum no.
9 4  no. of prim. and contr. sfunctions
4232.61  sexponents
634.882
146.097
42.4974
14.1892
1.9666
5.1477
0.4962
0.1533
.002029 .0 .0 .0  scontraction matrix
.015535 .0 .0 .0
.075411 .0 .0 .0
.257121 .0 .0 .0
.596555 .0 .0 .0
.242517 .0 .0 .0
.0 1.0 .0 .0
.0 .0 1.0 .0
.0 .0 .0 1.0
5 2  no. of prim. and contr. pfunctions
18.1557  pexponents
3.98640
1.14290
0.3594
0.1146
.018534 .0  pcontraction matrix
.115442 .0
.386206 .0
.640089 .0
.0 1.0
1 1  no. of prim. and contr. dfunctions
.75  dexponents
1.0  dcontraction matrix
C1 0.00000 0.00000 0.00000  atomlabel, Cartesian coordinates
C2 1.00000 0.00000 0.00000  atomlabel, Cartesian coordinates
End Of Basis  end of basis set definition
The label within the ECP library is given as input in the line following the keyword BASIS SET. The label defines either the valence basis set and core potential which is assigned to a frozencore atom or the embedding potential which is assigned to an environmental frozeion. Here, all the comments made about this label in the section The basis set label and the basis set library for allelectron basis sets stand, except for the following changes:
Various Density Functional Theory (DFT) models can be used in MOLCAS . Energies and analytical gradients are available for all DFT models. In DFT the exact exchange present in HF theory is replaced by a more general expression, the exchangecorrelation functional, which accounts for both the exchange energy, E_{X} [P] and the electron correlation energy ,E_{C} [P].
(6.15) 
The various DFT methods differ in which function, f, is used for E_{X}[P] and for E_{C}[P]. In MOLCAS pure DFT methods are supported, together with hybrid methods, in which the exchange functional is a linear combination of the HF exchange and a functional integral of the above form. The latter are evaluated by numerical quadrature. In the SEWARD input the parameters for the numerical integration can be set up. In the SCF and RASSCF inputs the keywords for using different functionals can be specified. Names for the various pure DFT models are given by combining the names for the exchange and correlation functionals.
The DFT gradients has been implemented for both the fixed and the moving grid approach [106,107,108]. The latter is known to be translationally invariant by definition and is recommended in geometry optimizations.
Below follows a description of the input to the numerical integration utility in the SEWARD input.
Compulsory keywords
Keyword  Meaning 
GRID Input  This marks the beginning of the input to the numerical integration utility. 
END Of GridInput  This marks the end of the input to the numerical integration utility. 
Optional keywords
Keyword  Meaning 
GRID  It specifies the quadrature quality. The possible indexes that can follow are COARSE, SG1GRID, FINE, ULTRAFINE following the Gaussian98 convention. Default is FINE. 
RQUAd  It specifies the radial quadrature scheme. Options are LOG3 (Mura and Knowles), BECKE (Becke) , MHL (Murray et a.), TA (Treutler and Ahlrichs, defined for HKr), and LMG (Lindh et al.), respectively. The default is MHL. 
GGL  It activates the use of Gauss and GaussLegendre angular quadrature. Default is to use the Lebedev angular grid. 
LEBEdev  It turns on the Lebedev angular grid. 
LOBAtto  It activates the use of Lobatto angular quadrature. Default is to use the Lebedev angular grid. 
LMAX  It specifies the angular grid size. Default is 29. 
NGRId  It specifies the maximum number of grid points to process at one instance. Default is 5500 grid points. 
NOPRunning  It turns off the the angular prunning. Default is to prune. 
NR  It is followed by the number of radial grid points. Default is 75 radial grid points. 
FIXEd grid  Use a fixed grid in the evaluation of the gradient. This corresponds to using the grid to numerically evaluate the analytic gradient expression. Default is to use a moving grid. 
MOVIng grid  Use a moving grid in the evaluation of the gradient. This correspond to evaluating the gradient of the numerical expression of the DFT energy. This is the default. 
THREshold  It is followed by the value for the the radial threshold. Default value is 1.0D13. 
T_X  Threshold for screening in the assembling of the density on the grid. Default value is 1.0D18. 
T_Y  Threshold for screening in the assembling of the integrals. Default value is 1.0D11. 
NOSCreening  Turn off any screening in the numerical integration. 
CROWding  The crowding factor, according to MHL, used in the pruning of the angular grid close to the nuclei. Default value 3.0. 
The SCF and RASSCF programs have their own keywords to decide which functionals to use in a DFT calculation.
Below follows an example of a DFT calculation with two different functionals.
&GATEWAY
Basis set
H.321G.....
H1 0.0 0.0 0.0
End of basis
&SEWARD
Grid input
RQuad= Log3; nGrid= 50000; GGL; lMax= 26; Global
End of Grid Input
&SCF; Occupations=1; KSDFT=LDA5; Iterations= 1 1
&SCF; Occupations= 1; KSDFT=B3LYP; Iterations= 1 1
The current different implementation of all relativistic operators in MOLCAS as described in the following subsubsections has been programmed and tested in Ref.[109]
For allelectron calculations, the preferred way is to use the scalarrelativistic DouglasKrollHess (DKH) Hamiltonian, which, in principle, is available up to arbitrary order in Molcas; for actual calculations, however, the standard 2nd order is usually fine, but one may use a higher order that 8th order by default to be on the safe side.
The arbitraryorder Hamiltonian is activated by setting
RXXPyy
somewhere in the SEWARD input, where the XX denotes the order of the DKH Hamiltonian in the external potential. I.e., for the standard 2ndorder Hamiltonian you may use R02O. Note in particular that the parametrization P does not affect the Hamiltonian up to fourth order. Therefore, as long as you run calculations with DKH Hamiltonians below 5th order you may use any symbol for the parametrization as they would all yield the same results.
The possible parametrizations P of the unitary transformation used in the DKH transformation supported by MOLCAS are:
Up to fourth order (XX=04) the DKH Hamiltonian is independent of the chosen parametrization. Higherorder DKH Hamiltonians depend slightly on the chosen parametrization of the unitary transformations applied in order to decouple the Dirac Hamiltonian. Since the EXP parameterization employs a fast algorithm [110], it is recommended for highorder DKH transformation.
For details on the arbitraryorder DKH Hamiltonians see [111] with respect to theory, [112] with respect to aspects of implementation, and [113] with respect to general principles of DKH. The current version of MOLCAS employs different algorithms, but the polynomial cost scheme of the DKH implementation as described in [110] is used as the default algorithm. The implementation in MOLCAS has been presented in [109].
For details on the different parametrizations of the unitary transformations see [114].
As mentioned above, fourcomponent molecular property operators need to be DKH transformed as well when going from a fourcomponent to a two or onecomponent description; the results do not coincide with the wellknown corresponding nonrelativistic expressions for a given property but are properly picture change corrected.
The transformation of electricfieldlike molecular property operators can be carried out for any order smaller or equal to the order chosen for the scalarrelativistic DKH Hamiltonian. In order to change the default transformation of order 2, you may concatenate the input for the DKH Hamiltonian by 2 more numbers specifying the order in the property,
RxxPyy
where yy denotes the order of the Hamiltonian starting with first order 01. The DKH transformation is then done automatically for all oneelectron electricfieldlike oneelectron property matrices.
Also note that the current implementation of both the Hamiltonian and the property operators is carried out in the full, completely decontracted basis set of the molecule under consideration. The local nature of the relativistic contributions is not yet exploited and hence large molecules may require considerable computing time for all higherorder DKH transformations.
For details on the arbitraryorder DKH properties see [115] with respect to theory and [116,109] with respect to implementation aspects.
Exact decoupling of the relativistic Dirac Hamiltonian can be achieved with infiniteorder approaches, such as the socalled X2C (exacttwocomponent) and BSS (BaryszSadlejSnijders) method. In Molcas, both methods are available for allelectron calculations. The evaluation of transformation matrices employs a noniterative scheme.
The exact decoupling Hamiltonian is activated by setting either RX2C or RBSS somewhere in the SEWARD input, where RX2C and RBSS denote using the scalar (onecomponent) X2C and BSS Hamiltonian respectively. The oneelectron Hamiltonian as well as the property integrals will be transformed according to the given exact decoupling method. In other words, all property integrals are by default picture change corrected.
The computation time of the X2C/BSS method is almost the same as of the DKH method at 8th order, while X2C is a little bit faster than BSS since the additional freeparticle FoldyWouthuysen transformation is skipped in the X2C approach[109]. For molecules including only light atoms, the DKH method with low orders (<8) is enough to account for the relativistic effects.
The differences between different exact decoupling relativistic methods are very small compared to errors introduced by other approximations, such as the basis set incompleteness, approximate density functionals, etc. Therefore, any exact decoupling model is acceptable for the treatment of relativistic effects in molecular calculations.
For details on the exact decoupling approaches see [109] with respect to theories and comparison of numerical results, [117,118,119] for the X2C method, and [120,121] for the BSS method.