PAObased local correlation treatments
Introduction
In this section the original PAObased local correlation methods are described. Since Molpro version 2018, a completely new PNOLCCSD(T)F12 program is available, which is much more accurate and well parallelized. It is recommended to use this new program, which is described in section local correlation methods with pair natural orbitals (PNOs).
The PAO based local correlation program of Molpro can currently perform closedshell LMP2, LMP3, LMP4(SDTQ), LCISD, LDCSD, LQCISD(T), and LCCSD(T) calculations. For large molecules, all methods scale linearly with molecular size, provided very distant pairs are neglected, and the integraldirect algorithms are used.
Much higher efficiency is achieved by using density fitting (DF) approximations to compute the integrals. Density fitting is available for all local methods up to LCCSD(T), as well as for analytical LMP2 gradients. Only iterative triples methods like LCCSDT1b can currently not be done with density fitting.
The errors introduced by DF are negligible, and the use of the DF methods is highly recommended. Linear scaling can be obtained in DFLMP2 using the LOCFIT
option (see Ref. 11); in DFLCCSD(T), the most important parts also scale linearly, but some transformation steps scale quadratically.
Energy gradients are available for LMP2, DFLMP2, DFSCSLMP2, and LQCISD (in the latter case only for LOCAL=1
, i.e. the local calculation is simulated using the canonical program, and savings only result from the reduced number of pairs).
Local explicitly correlated methods (DFLMP2R12 and DFLMP2F12, DFLCCSD(T)F12 are described in section explicitly correlated methods.
Before using these methods, it is strongly recommended to read the literature in order to understand the basic concepts and approximations. A recent review [1] and Ref. [2] may be suitable for an introduction.
References:
Review:
$[1]$ H.J. Werner and K. Pflüger, On the selection of domains and orbital pairs in local correlation treatments, Ann. Rep. Comp. Chem. 2, 53 (2006).
General local Coupled Cluster:
$[2]$ C. Hampel and H.J. Werner, Local Treatment of electron correlation in coupled cluster (CCSD) theory, J. Chem. Phys. 104, 6286 (1996).
$[3]$ M. Schütz and H.J. Werner, Local perturbative triples correction (T) with linear cost scaling, Chem. Phys. Lett. 318, 370 (2000).
$[4]$ M. Schütz, Loworder scaling local electron correlation methods. III. Linear scaling local perturbative triples correction (T), J. Chem. Phys. 113, 9986 (2000).
$[5]$ M. Schütz and H.J. Werner, Loworder scaling local electron correlation methods. IV. Linear scaling local coupledcluster (LCCSD), J. Chem. Phys. 114, 661 (2001).
$[6]$ M. Schütz, Loworder scaling local electron correlation methods. V. Connected Triples beyond (T): Linear scaling local CCSDT1b, J. Chem. Phys. 116, 8772 (2002).
$[7]$ M. Schütz, A new, fast, semidirect implementation of Linear Scaling Local Coupled Cluster Theory, Phys. Chem. Chem. Phys. 4, 3941 (2002).
Multipole treatment of distant pairs:
$[8]$ G. Hetzer, P. Pulay, H.J. Werner, Multipole approximation of distant pair energies in local MP2 calculations, Chem. Phys. Lett. 290, 143 (1998).
Linear scaling local MP2:
$[9]$ M. Schütz, G. Hetzer and H.J. Werner, Loworder scaling local electron correlation methods. I. Linear scaling local MP2, J. Chem. Phys. 111, 5691 (1999).
$[10]$ G. Hetzer, M. Schütz, H. Stoll, and H.J. Werner, Loworder scaling local electron correlation methods II: Splitting the Coulomb operator in linear scaling local MP2, J. Chem. Phys. 113, 9443 (2000).
Density fitted local methods:
$[11]$ H.J. Werner, F. R. Manby, and P. J. Knowles, Fast linear scaling secondorder MøllerPlesset perturbation theory (MP2) using local and density fitting approximations, J. Chem. Phys. 118, 8149 (2003). $[12]$ M. Schütz and F.R. Manby, Linear scaling local coupled cluster theory with density fitting. Part I: 4 external integrals, Phys. Chem. Chem. Phys. 5, 3349 (2003).
$[13]$ Polly, H.J. Werner, F. R. Manby, and Peter J. Knowles, Fast HartreeFock theory using local density fitting approximations, Mol. Phys. 102, 2311 (2004).
$[14]$ H.J. Werner and M. Schütz, Loworder scaling coupled cluster methods (LCCSD(T)) with local density fitting approximations, in preparation.
LMP2 Gradients and geometry optimization:
$[15]$ A. El Azhary, G. Rauhut, P. Pulay and H.J. Werner, Analytical energy gradients for local secondorder MøllerPlesset perturbation theory, J. Chem. Phys. 108, 5185 (1998).
$[16]$ G. Rauhut and H.J. Werner, Analytical Energy Gradients for Local CoupledCluster Methods, Phys. Chem. Chem. Phys. 3, 4853 (2001).
$[17]$ M. Schütz, H.J. Werner, R. Lindh and F.R. Manby, Analytical energy gradients for local secondorder MøllerPlesset perturbation theory using density fitting approximations, J. Chem. Phys. 121, 737 (2004).
LMP2 vibrational frequencies:
$[18]$ G. Rauhut, A. El Azhary, F. Eckert, U. Schumann and H.J. Werner, Impact of Local Approximations on MP2 Vibrational Frequencies, Spectrochimica Acta 55, 651 (1999).
$[19]$ G. Rauhut and H.J. Werner The vibrational spectra of furoxan and dichlorofuroxan: a comparative theoretical study using density functional theory and Local Electron Correlation Methods, Phys. Chem. Chem. Phys. 5, 2001 (2003).
$[20]$ T. Hrenar, G. Rauhut and H.J. Werner, Impact of local and density fitting approximations on harmonic vibrational frequencies, J. Phys. Chem. A., 110, 2060 (2006).
Intermolecular interactions and the BSSE problem:
$[21]$ M. Schütz, G. Rauhut and H.J. Werner, Local Treatment of Electron Correlation in Molecular Clusters: Structures and Stabilities of (H$_2$O)$_n$, $n=24$, J. Phys. Chem. 102, 5997 (1998). See also [2] and references therein.
$[22]$ N. Runeberg, M. Schütz and H.J. Werner, The aurophilic attraction as interpreted by local correlation methods, J. Chem. Phys. 110, 7210 (1999).
$[23]$ L. Magnko, M. Schweizer, G. Rauhut, M. Schütz, H. Stoll, and H.J. Werner, A Comparison of the metallophilic attraction in (XMPH$_3$)$_2$ (M=Cu, Ag, Au; X=H, Cl), Phys. Chem. Chem. Phys. 4, 1006 (2002).
Improved treatment of intermolecular pairs in local coupled cluster methods (beyond LMP2):
$[24]$ O. Masur, D. Usvyat and M. Schütz Efficient and accurate treatment of weak pairs in local CCSD(T) calculations, J. Chem. Phys. 139, 164116 (2013).
$[25]$ M. Schütz, O. Masur and D. Usvyat Efficient and accurate treatment of weak pairs in local CCSD(T) calculations. II. Beyond the ring approximation, J. Chem. Phys. 140, 244107 (2014).
Getting started
The local correlation treatment is switched on by preceding the command name by an L
, i.e., by using the LMP2
, LMP3
, LMP4
, LQCISD
, LCCSD
, LDCSD
, or LCISD
commands.
The LQCISD
and LCCSD
commands can be appended by a specification for the perturbative treatment of triple excitations (e.g., LCCSD(T0)
):
(T)
Use the default triples method. Currently this is T0.(T0)
Noniterative local triples. This is the fastest triples option. It is usually sufficiently accurate and recommended to be used in most cases.(T1)
T0 plus one perturbative update of the triples amplitudes. If the accuracy of T0 is insufficient (very rarely the case!), this can be used to improve the accuracy. The computational cost is at least twice as large as for T0. In contrast to T0, the triples amplitudes must be stored on disk, which can be a bottleneck in calculations for large molecules. Also the memory requirements are substantially larger than for T0.(T1C)
As T1, but a caching algorithm is used which avoids the simultaneous storage of all triples amplitudes on disk (as is the case for(T1)
or(TF)
). Hence, T1C requires less disk space but more CPUtime than T1. The more disk space is made available for the caching algorithm (using the T1DISK option on the local card, see below), the less CPU time is used.(TF)
Full iterative triples calculation. With full domains and without weak pair approximations this gives the same result as a canonical (T) calculation. Typically, 35 iterations are needed, and therefore the computational effort is 23 times larger than for (T1). The disk and memory requirements are the same as for T1. The T0 energy is also computed and printed.TFULL
andFULL
are aliases forTF
.(TA)
As TF, but the T1 energy is also computed. Since the first iteration is different for T1, the convergence of the triples iterations is slightly different with TF and TA (TF being somewhat faster in most cases).TALL
andALL
are aliases forTA
.
Density fitting can be invoked by prepending the command name by DF
, e.g. DFLMP2
, DFLCCSD(T0)
etc. In density fitting calculations an additional auxiliary basis set is needed. Details about choosing such basis sets and other options for density fitting are described in sections Densityfitted LMP2 (DFLMP2), LDCSD (DFLDCSD), and coupled cluster (DFLCCSD(T0)) and density fitting.
The general input for local coupled LMP2 or coupled cluster calculations is:
LMP2
,options Local MP2 calculationLCCSD
,options Local CCSD calculationLDCSD
,options Local DCSD calculationLCCSD(T0)
,options Local CCSD(T0) calculation
The same options as on the command line can also be given on subsequent LOCAL
and MULTP
directives. Instead of using the MULTP
directive, the MULTP
option on the command line can also be used.
In the following, we will first give a summary of all options and directives. These will be described in more detail in the subsequent sections. For new users it is recommended to read section doing it right at the end of this chapter before starting calculations.
Summary of options
Many options can be specified on the command line. For all options appropriate default values are set, and so these options must usually be modified only for special purposes. For convenience and historical reasons, alias names are available for various options, which often correspond to the variable name used in the program. The following table summarizes the options, aliases and default values. In the following, the parameters will be described in more detail.
Summary of local (multp ) options and their default values 


Parameter  Alias  Default value  Meaning 
General Parameters  
LOCAL  4  determines which program to use.  
MULTP  0  turns on multipole approximations for distant pairs.  
SAVEDOM  SAVE  0  specifies record for saving domain info. 
RESTDOM  START  0  specifies record for reading domain info. 
LOCORB  0  activates or deactivates orbital localization.  
LOC_METHOD  specifies which localization method to use.  
CANONICAL  0  allows to use canonical virtual orbitals (for testing).  
PMDEL  CPLDEL  0  discards contributions of diffuse functions in PM localization. 
SAVORB  SAVLOC  0  specifies record for saving local orbitals. 
DOMONLY  0  if 1, only domains are made. if 2, only orbital domains are made.  
Parameters to define domains  
THRBP  DOMSEL  0.98  BoughtonPulay selection criterion for orbital domains. 
NPASEL  0  charge used in the NPA selection criterion for orbital domains.  
CHGMIN  0.01  determines the minimum allowed atomic charge in domains.  
CHGMINH  0.03  as CHGMIN , but used for Hatoms (default 0.03). 

CHGMAX  0.40  If the atomic charge is larger than this value, the atom is always included in the domain.  
MAXANG  MAXL  99  angular momentum restriction for BP domain selection 
MAXBP  0  determines how atoms are ranked in BP procedure.  
MULLIKEN  LOCMUL  0  determines the method to determine atomic charges. 
MERGEDOM  0  merges overlapping domains.  
DELCOR  IDLCOR  2  delete projected core AOs up to certain shell. 
DELBAS  IBASO  0  determines how to remove redundancies. 
Distance criteria for domain extensions  
REXT  0  criterion for all pair domains.  
REXTS  0  criterion for strong pair domains.  
REXTC  0  criterion for strong and close pair domains.  
REXTW  0  criterion for strong, close, and weak pair domains.  
Connectivity criteria for domain extensions  
IEXT  0  criterion for all pair domains.  
IEXTS  0  criterion for strong pair domains.  
IEXTC  0  criterion for strong and close pair domains.  
IEXTW  0  criterion for strong, close, and weak pair domains.  
Parameters to select pair classes  
USE_DIST  1  determines if distance of connectivity criteria are used.  
RCLOSE  CLOSEP  1  distance criterion for selection of weak pairs. 
RWEAK  WEAKP  3  distance criterion for selection of weak pairs. 
RDIST  DISTP  8  distance criterion for selection of distant pairs. 
RVDIST  VERYD  15  distance criterion for selection of very distant pairs. 
ICLOSE  1  connectivity criterion for selection of weak pairs.  
IWEAK  2  connectivity criterion for selection of weak pairs.  
IDIST  5  connectivity criterion for selection of distant pairs.  
IVDIST  8  connectivity criterion for selection of very distant pairs.  
CHGMIN_PAIRS  CHGMINP  0.20  determines minimum charge of atoms used for pair classification. 
KEEPCL  0  determines if close pairs are included in LCCSD.  
Parameter for multipole treatment of exchange operators  
DSTMLT  3  multipole expansion level for distant pairs  
Parameters for energy partitioning  
IEPART  0  If nonzero: do energy partitioning.  
EPART  3.0  cutoff parameter for determining individual monomers.  
Parameters for redundancy check using DELBAS=1 (not recommended) 

TYPECHECK  TYPECHK  1  activates basis function type restrictions. 
DELSHL  IDLSHL  1  determines if whole shells are to be deleted. 
DELEIG  IDLEIG  1  determines how to select redundant functions. 
DELCMIN  CDELMIN  0.1  parameter for use with DELEIG=1 
Parameters for choosing operator domains in LCCSD  
OPDOM  IOPDOM  5  determines how operator domains are determined for LCCSD 
RMAXJ  8  distance criterion for Joperator list.  
RMAXK  8  distance criterion for Koperator list.  
RMAXL  15  distance criterion for Loperator list.  
RMAX3X  5  distance criterion for 3ext integral list.  
RDOMJ  0  distance criterion for Koperator domains.  
RDOMK  8  distance criterion for Joperator domains.  
IMAXJ  5  connectivity criterion for Joperator list.  
IMAXK  5  connectivity criterion for Koperator list.  
IMAXL  8  connectivity criterion for Loperator list.  
IMAX3X  3  connectivity criterion for 3ext integral list.  
IDOMJ  0  connectivity criterion for Koperator domains.  
IDOMK  5  connectivity criterion for Joperator domains.  
Miscellaneous options  
SKIPDIST  SKIPD  3  determines at which stage weak and distant pairs are eliminated 
ASYDOM  JITERM  0  parameter for use of asymmetric domains 
LOCSING  LOCSNG  0  determines virtual space used for singles 
PIPEKAO  LOCAO  0  activates AO localization criterion 
NONORM  2  determines whether projected functions are normalized  
LMP2ALGO  MP2ALGO  3  if nonzero, use loworder scaling method in LMP22 iterations 
OLDDEF  0  allows to revert to older defaults  
T1DISK  10  maximum disk space (in GByte) for T1 caching algorithm  
Thresholds  
THRBP  0.98  Threshold BoughtonPulay method.  
THRPIP  1.d9  Threshold for PipekMezey or Boys localization.  
THRCPL  1.d9  Threshold for coupledperturbed localization  
THRORB  1.d6  Threshold for eliminating projected orbitals with small norm.  
THRLOC  1.d6  Threshold for eliminating redundant projected orbitals.  
THRCOR  1.d1  Threshold for eliminating projected core orbitals.  
THRMP2  1.d8  Threshold for neglecting small fock matrix elements in the LMP2 iteration. 
Summary of directives
The same standard directives as in the canonical programs, e.g., OCC
, CLOSED
, CORE
, WF
, ORBITAL
are also valid in the local methods. In addition, there are some directives which only apply to local calculations:
LOCAL
Invokes local methods and allows to specify the same options as on the command line.MULTP
AsLOCAL
, but multipole approximations are used for distant pairs.DOMAIN
Define domains manually (not recommended).MERGEDOM
Allows to merge domainsREGION
Allows to select regions of a molecule to be treated at a certain level of theory.ENEPART
Analysis of pair energies.SAVE
Save domains and LCCSD amplitudes.START
Restart with domains and LCCSD amplitudes from a previous calculation.
General Options
LOCAL
=local Determines which method is used:
LOCAL=0
: Conventional (nonlocal) calculation.
LOCAL=1
: Local method is simulated using canonical MOs. The local basis is used only at an intermediate stage to update the amplitudes in each iteration (only for testing).
LOCAL=2
: Calculation is done in local basis, but without using local blocking (i.e. full matrices are used). This is the most expensive method and only for testing.
LOCAL=3
: Fully local calculation (obsolete).
LOCAL=4
: Fully local calculation (default). This is the fastest method for large molecules with many weak pairs and requires minimum memory.
LOCORB
=option If this option is given and option$\gt 0$, the orbitals are localized using the PipekMezey technique. If this option is not given or option=0 (default), the orbitals are localized unless localized orbitals are found in the orbital record (cf.ORBITAL
directive andLOCALI
command). In the latter case, the most recent localized orbitals are used. Setting option=1 switches the localization off. If option$\gt 1$ the localized orbitals are printed. Note: Boys localization can only be performed using theLOCALI
command. The program will use the Boys orbitals if they are found in the orbital record and theLOCORB
option is absent or option$\le 0$.LOC_METHOD
=method This option allows to select between PipekMezey (method=PM
) or Natural Orbitals (method=NBO
) localization. If PipekMezey orbitals are requested, the default BoughtonPulay domain selection will be used. When method=NBO
, the domain selection will be based on the NPA charges, withNPASEL
=0.03 (by default). In both cases, the domain selection parameter can be explicitly given (cf.DOMSEL
andNPASEL
options) in order to use another domain criterion. If localized orbitals are found in the orbital record, but the type does not coincide, the orbitals are again localized according to method and used in the calculation.SAVORB
=record Allows the localized and projected orbitals to be saved in record=name.ifil for later use (e.g. plotting). The two orbital sets are stored in the same dump record and can be restored at later stages using
ORBITAL
,record,[TYPE
=]LOCAL
or
ORBITAL
,record,[TYPE
=]PROJECTED
, respectively.
DOMONLY
=value If value$\gt 0$ only domains are made, but no energy is computed. This can be used to check and save the domains for later use.DSTMLT
=level Determines the expansion level of the multipole expansion of distant pairs (e.g. 1 means dipole approximation, 2 quadrupole approximation and so on). The default forMULTP
is 3.INTERACT
Automatically determine individual molecules in a calculation and set appropriate pair lists for computing interaction energies. See section intermolecular interactions for more details. Parameters for energy partitioning:
IEPART
=value enables/disables energy partitioning.
iepart=0: Energy partitioning is disabled.
iepart=1: Energy partitioning is enabled.
iepart=2: Energy partitioning is enabled. Additionally, a list of all pair energies and their components is printed.
EPART
=cutoff Cutoff parameter to determine individual monomers in a cluster (i.e. centre groups). Should be somewhat larger than the largest intramolecular bond length (given in a.u.).
 Miscellaneous options:
SKIPDIST
=skipdist Testparameter. Its value should only affect the efficiency but not influence the results.
skipdist=1: Weak and distant pairs are set to zero after MP2 but are not eliminated from the pair list and not skipped in any loop.
skipdist=0: No pairs are deleted from pair list, but weak and distant pairs are skipped in the loops were appropriate.
skipdist=1: Very distant pairs are neglected from the beginning. Distant pairs are eliminated after MP2.
skipdist=2: As skipdist=1, but also weak pairs are eliminated after MP2.
skipdist=3: As skipdist=2, but distant pairs are eliminated from the operator list in case of LMP2 with multipole approximations for distant pairs. This is the default.
ASYDOM
=jiterm Experimental test parameter. Enables the use of asymmetric domains for distant pairs. The asymmetric domain approximation supplements the multipole approximation for distant pairs, as it suppresses the treatment of configurations for which no integrals can be computed by multipole expansion. This leads to computational savings and improved numerical stability.
jiterm=0: Disable asymmetric domains.
jiterm=1: Enable asymmetric domains (default).
jiterm=2: Enable a variation of the asymmetric domain formalism: Exchange operators will initially be projected to the asymmetric domain instead of simply packed.
LOCSING
=locsing If locsing.ne.0, the single excitations use the full space, i.e., they are not treated locally. This is only works forLOCAL=1
.MAXANG
=lmax The purpose of this experimental option is to reduce the basis set sensitivity of the BoughtonPulay (BP) method for domain selection. Only basis functions with angular momentum up to lmax1 are included when computing the overlap of the approximate and exact orbitals. For example,MAXANG=2
means to omit all contributions of $d$, $f$ and higher angular momentum functions. To obtain reasonable domains, the value ofTHRBP
must often be reduced (to 0.97 or so). This option should only be used with care!PIPEKAO
=option If option$\ge 0$, the orbitals are localized my maximizing the coefficients of basis functions of a given type at a given atom. Normally, this is only useful to uniquely define degenerate orbitals in atoms. For instance, when this option is used to localize the orbitals for a dimer like (Ar)$_2$ at a very long distance, clean $s$, $p_x$, $p_y$, and $p_z$ atomic orbitals will be obtained. It is not recommended to use this option for molecular calculations!NONORM
=value Determines if projected functions are normalized (not recommended). value=1: projected orbitals are normalized before redundancy check.
value=0: projected orbitals are normalized after redundancy check (default).
value=1: projected orbitals are normalized in redundancy check, afterwards unnormalized.
value=2: projected orbitals are never normalized (default in gradient calculations).
LMP2ALGO
=value If nonzero, use loworder scaling method in LMP2 iterations. Values can be 1, 2, or 3, and 3 is usually fastest if large basis sets are used.OLDDEF
=value For compatibility with older versions: if nonzero, revert to old defaults. Options set before this may be overwritten. Thresholds:
THRPIP
=thresh Threshold for PipekMezey or Boys localization. The localization is assumed to be converged if all $2 \times 2$ rotation angles are smaller than thresh. The default is $1.d9$. It can also be modified globally usingGTHRESH, LOCALI=
thresh.THRCPL
=thresh Threshold for coupled perturbed localization (default 1.d9). This is only needed in local gradient calculations.THRORB
=thresh Threshold for eliminating functions from pair domains whose norm is smaller than thresh after projecting out the occupied space. The default is throrb=1.d6.THRLOC
=thresh Threshold for eliminating redundant basis functions from pair domains. For each eigenvalue of ${\bf \tilde S}^{ij} \lt thresh$ one function is deleted. The default is 1.d6. The method used for deleting functions depends on the parametersIDLEIG
andIBASO
.THRMP2
=thresh Threshold for neglecting small fock matrix couplings in the LMP2 iterations (default 1.d8). Specifying a larger threshold speeds up the iterations but may lead to small errors in the energy. In the initial iterations, a larger threshold is chosen automatically. It is gradually reduced to the specified final value during the iterations.THRCOR
=thresh Threshold for deleting projected core orbitals. The functions are only deleted if their norm is smaller than thresh (default 0.1)
The thresholds can also be specified on the THRESH
directive.
Options for selection of domains
The following sections describe the most important options which affect the domains.
Standard domains
Standard domains are always determined first. They are used to define strong, close, weak, and distant pairs. More accurate results can be obtained with extended domains, as described in section extended domains.
THRBP
=value Threshold for selecting the atoms contributing to orbital domains using the method of Boughton and Pulay (BP). As many atoms as needed to fulfill the BP criterion are included in a domain. The order in which atoms are considered depends on the parameterMAXBP
, see below. The default isTHRBP=0.98
.THRBP=1.0
includes all atoms into each orbital domain, i.e., leads to full domains. If no pairs are neglected, this should yield the canonical MP2 energy. The criterion is somewhat basis dependent. See section orbital domains for recommended values of this threshold.CHGMIN
=value determines the minimum allowed Mulliken (or Löwdin) charge for an atom (except H) in a domain, i.e., atoms with a smaller (absolute) charge are not included, even if theTHRBP
criterion is not fulfilled (default 0.01).CHGMINH
=value asCHGMIN
, but used for Hatoms (default 0.03).CHGMAX
=value If the atomic charge is larger than this value, the atom is included, independent of any ranking.MAXBP
=maxbp If maxbp=1, the atoms are ranked according to their contribution to the BoughtonPulay overlap. If maxbp=0 (default), the atoms are ranked according to atomic charges. In both cases atoms with charges greater thanCHGMAX
are always included, and atoms with the same charges are added as groups.MULLIKEN
=option Determines the method to determine atomic charges.MULLIKEN=0
(default): squares of diagonal elements of ${\bf S}^{\frac{1}{2}} {\bf C}$ are used (Löwdin charges);MULLIKEN=1
: Mulliken gross charges are used. The first choice is less basis set dependent and more reliable with diffuse basis sets.MERGEDOM
=number If number is greater than zero, all orbital domains containing number or more atoms in common are merged (number=1 is treated as number=2, default 0). This is particularly useful for geometry optimizations of conjugated or aromatic systems like, e.g., benzene. In the latter case,MERGEDOM=1
causes the generation of full $\pi$domains, i.e., the domains for all three $\pi$orbitals comprise all carbon basis functions. Note that the merged domains are generated after the above print of orbital domains, and information about merged domains is printed separately. See section gradients and frequency calculations for further discussion of geometry optimizations.
There are some other options which should normally not be modified:
This parameter determines the method for eliminating redundant functions of pair domains.
ibaso=0: The space of normalized eigenvectors of ${\bf \tilde S}^{ij}$, which correspond to small eigenvalues, is eliminated (default). Any other value is not recommended and not further documented.
Activates elimination of basis functions corresponding to core orbitals. If nshell=1, only $1s$functions are eliminated from projected space. If nshell=2 (default) $1s$ functions on firstrow atoms, and $1s$, $2s$, and $2p$functions are eliminated on secondrow atoms. Nothing is eliminated on H or He atoms. If effective core potentials are used, nothing is deleted at the corresponding atom. Also, functions are only deleted if the norm of the projected function is below THRCOR
(default 0.1)
Extended domains
There are two alternative modes for domain extensions: either distance criteria REXT
, REXTS
, REXTC
, or REXTW
can be used. These are in Bohr and refer to the minimum distance between any atom in a standard orbital domain $[ij]$ and another atom. If an atom is found within the given distance, all PAOs at this atom are added to the domain $[ij]$. Alternatively, connectivity criteria IEXT
, IEXTS
, IEXTC
, or IEXTW
can be used. These refer to the number of bonds between any atom contained in the standard domain $[ij]$ and another atom. The advantage of distance criteria is that they select also atoms within the given radius which are not connected to the present domain by bonds. On the other hand, the connectivity criteria are independent of different bond lengths, e.g., for first and secondrow atoms. Only one of the two possibilities can be used, i.e., they are mutually exclusive.
REXT
=value Distance criterion for extension of all pair domains.REXTS
=value Distance criterion for extension of strong pair domains.REXTC
=value Distance criterion for extension of strong and close pair domains.REXTW
=value Distance criterion for extension of strong, close, and weak pair domains.IEXT
=value Connectivity criterion for extension of all pair domains.IEXTS
=value Connectivity criterion for extension of strong pair domains.IEXTC
=value Connectivity criterion for extension of strong and close pair domains.IEXTW
=value Connectivity criterion for extension of strong, close, and weak pair domains.
By default, domains are not extended, i.e., the default values of all parameters listed above are zero. Note that the pair classes are determined on the basis of the standard domains, and therefore domain extensions have no effect on the pair lists.
Also note that the computational effort increases with the fourth power of the domain sizes and can therefore increase quite dramatically when extending domains. This does not affect the linear scaling behaviour in the asymptotic limit.
Manually Defining orbital domains (DOMAIN)
It is possible to define the domains “by hand”, using the DOMAIN
directive:
DOMAIN
,orbital,atom1, atom2 …
where orbital has the form iorb.isym
, e.g., 3.1 for the third orbital in symmetry 1, and atomi are the atomic labels as given in the Zmatrix geometry input, or, alternatively, the Zmatrix row numbers. All basis functions centred at the given atoms are included into the domain. For instance
DOMAIN,3.1,C1,C2
defines a domain for a bicentric bond between the carbon atoms C1
and C2
. The DOMAIN
directive must be given after any OCC
, CLOSED
, or CORE
directives. Note that the order of the localized orbitals depends on the localization procedure, and could even change as function of geometry, and therefore manual DOMAIN
input should be used with great care. The domains of all orbitals which are not explicitly defined using DOMAIN
directive are determined automatically as usual.
Options for selection of pair classes
There are two alternative modes for defining the pair classes: either distance criteria RCLOSE
, RWEAK
, RDIST
, RVDIST
can be used. These are in Bohr and refer for a given orbital pair $(ij)$ to the minimum distance $R^{(ij)}$ between any atom in the standard orbital domains $[i]$ and any atom in the standard orbital domains $[j]$ . Alternatively, the connectivity criteria ICLOSE
, IWEAK
, IDIST
, IVDIST
can be used. These refer to the minimum number of bonds between any atom contained in the standard domain $[i]$ and any atom contained in the standard domain $[j]$ The advantage of using connectivity criteria is the independence of the bond lengths, while the advantage of distance criteria (default) is that they are also effective in nonbonding situations. Only one of the two possibilities can be used, i.e., they are mutually exclusive. The use of distance criteria is the default. Using connectivity criteria for pair selection requires to set the option USE_DIST=0
.
USE_DIST
(default 1) If nonzero, use distance criteria, otherwise connectivity criteria.CHGMIN_PAIRS
Only atoms in the primary domains are considered for the pair classification if the atomic Löwdin charge is larger thanCHGMIN_PAIRS
(default value 0.2). This criterion was introduced in order to reduce the dependence of the pair selection on localization tails.RCLOSE
(default 1) Strong pairs are defined by $0 \le R^{(ij)} \lt {\tt RCLOSE}$. Close pairs are defined by ${\tt RCLOSE} \le R^{(ij)} \lt {\tt RWEAK}$.RWEAK
(default 3) Weak pairs are defined by ${\tt RWEAK} \le R^{(ij)} \lt {\tt RDIST}$.RDIST
(default 8) Distant pairs are defined by ${\tt RDIST} \le R^{(ij)} \lt {\tt RVDIST}$.RVDIST
(default 15) Very distant pairs for which $R^{(ij)} \ge$RVDIST
are neglected.ICLOSE
(default 1) Strong pairs are separated by less thatICLOSE
bonds. Close orbital pairs are separated by at leastICLOSE
bonds but less thanIWEAK
bonds.IWEAK
(default 2) Weak orbital pairs are separated by at leastIWEAK
bonds but less thanIDIST
bonds.IDIST
(default 5) Distant orbital pairs are separated by at leastIDIST
bonds but less thanIVDIST
bonds.IVDIST
(default 8) Very distant orbital pairs (neglected) are separated by at leastIVDIST
bonds.KEEPCL
(default 0) IfKEEPCL=1
, the LMP2 amplitudes of close pairs are included in the computation of the strong pair LCCSD residuals. IfKEEPCL=2
all close pairs are fully included in the LCCSD (this does not affect the triples list). This option is not yet implemented as efficiently as it could, and can therefore lead to a significant increase of the CPU time.
Setting RCLOSE
or RWEAK
to zero means that all pairs up to the corresponding class are treated as strong pairs (RWEAK=0
implies RCLOSE=0
). For instance, RCLOSE
=0 means that strong and close pairs are fully included in the LCCSD (in this case KEEPCL=1
has no effect). Note, however, that setting RCLOSE
=0 increases the length of the triples list. Setting RDIST=0
means that all distant pairs are treated as weak pairs. This does not affect RWEAK
and RCLOSE
and has no effect unless multipole approximations are used for distant pairs. Setting RVDIST=0
means that no very distant pairs are neglected. Again, this has no effect on the other distance parameters.
Directives
The LOCAL directive
The LOCAL
directive can be used to specify options for local calculations. If this directive is inside the command block of a local calculation, the options are used only for the current calculation, and this is entirely equivalent as if they were specified on the command line. The LOCAL
directive can also be given outside a command block, and in this case the options are used for all subsequent local correlation calculations in the same input.
Example:
DFLMP2,THRBP=0.985
is equivalent to
{DFLMP2 LOCAL,THRBP=0.985}
In the following example the LOCAL
directive is global and acts on all subsequent local calculations, i.e. both calculations will use THRBP=0.985
LOCAL,THRBP=0.985 DFLMP2 !local MP2 calculation OPTG !geometry optimization using the DFLMP2 energy DFLCCSD(T) !local coupled cluster at the optimized structure.
The MULTP directive
The MULTP
directive turns on the multipole approximations for distant pairs, as described in Ref. [8]. Further options can be given as described above for the LOCAL
directive.
LOCAL,MULTP
,options
is equivalent to
MULTP
,options
The level of the multipole approximation can be chosen using option DSTMLT
(default 3) ( 1 means dipole approximation, 2 quadrupole approximation and so on).
The multipole approximation reduces the computational cost of LMP2 calculations for very large molecules, but leads to some additional errors, see Ref. [8]. It is normally not recommended to be used in coupledcluster calculations and should never be used for computing intermolecular forces. It can also not be used in geometry optimizations or gradient calculations.
Saving the wavefunction (SAVE)
The wavefunction can be saved for later restart using
SAVE
,record
where record has the usual form, e.g., 4000.2 means record 4000 on file 2. If this directive is given, the domain information as well as the amplitudes are saved (for MPn the amplitudes are not saved). If just the domain information should be stored, the SAVE
option on the LOCAL
directive must be used (cf. section summary of options).
Restarting a calculation (START)
Local CCSD, DCSD or QCISD calculations can be restarted using
START
,record [sec:locstart]
The record given must have been saved in a previous local calculation using the SAVE
directive (otherwise this directive is ignored). If the START
directive is given, the domain information as well as the amplitudes of the previous calculation are used for restart. It is possible, for instance, to start a local CCSD calculation with the amplitudes previously saved for a local QCISD calculation (but of course it is not possible to use a record saved for a nonlocal CCSD or QCISD calculation). If it is intended only to use the domain information but not the amplitudes for a restart, the START
option on the command line or LOCAL
directive must be used (cf. section summary of options).
Correlating subsets of electrons (REGION)
In large molecules, it may be sufficient to correlate only the electrons in the vicinity of an active group, and to treat the rest of the molecule at the SCF level. This approach can even be extended, different correlation levels may be used for different sections of the system. The REGION
directive allows the specification of a subset of atoms:
REGION
,METHOD
=method,[DEFAULT
=default_method], [TYPE
=INCLUSIVEEXCLUSIVE
], atom1, atom2 …
The orbitals located at these atoms will be treated at the level specified in method. The remaining orbitals will be treated as defined in default. If not given by the user, the latter option will be set to HF
.
The orbital selection can be done in two ways. If type is set to INCLUSIVE
, any orbital containing one of the atoms in its domain centre list will be included. If type is set to EXCLUSIVE
, the program will only add orbitals whose domains are exclusively covered by the given atoms. Any local correlation treatment can be given as method, with the restriction that only MP2
and HF
can be used as default_method. Up to two REGION
directives may be included in a single calculation, ordered according to the correlation level (method) specified for the region. The highest level region should be given last.
It is advisable to check the region orbital list and the orbital domains printed by the program. The use of regions may significantly reduce the computation time, and, provided the active atoms are sensibly chosen, may give still sufficiently accurate results for the active group, e.g. bond lengths and bond angles.
Domain Merging (MERGEDOM)
The restriction of the virtual space in local calculations may result in discontinuities for reaction path calculations due to changes of the geometry dependent domains. This may be avoided by the use of a MERGEDOM
directive
MERGEDOM
,[NEIGHBOUR
=value],[CENTERS
=[atom1, atom2…]], [RECORD
=…],CHECK
This directive provides augmented domains, which can be saved (using option or directive SAVE
, see section saving the wavefunction (SAVE)) for later use in reaction paths or in single point calculations (in cases where the orbital domain description is unbalanced). The use of the neighbour option works in the same way as the local option MERGEDOM
, with value specifying the number of coincident centres. If the centres option is used, an atom list should be given (enclosed by square brackets). The domains of all orbitals located exclusively at these atoms will be merged, and the resulting merged domains will be used for all these orbitals.
One may also give a record number from a previously saved local calculation. The domain list contained in the record will be matched to the current one, and orbital domains augmented (merged) to include both sets. This domain definition should then be adequate for calculations on both points (and all those in between). This procedure can be repeated to include more geometries. In this way domains can be defined that are appropriate for a whole range of geometries (e.g. a reaction path), and if these domains are used in all calculations a strictly smooth potential energy surface is obtained.
Energy partitioning for molecular cluster calculations (ENEPART)
The local character of occupied and virtual orbitals in the local correlation treatment also offers the appealing possibility to decompose the intermolecular interaction energy of molecular clusters into individual contributions of different excitation classes. This allows to distinguish between intramolecular, dispersive, and ionic components of the correlation contribution to the interaction energy (cf. M. Schütz, G. Rauhut and H.J. Werner, J. Phys. Chem. 102, 5197 (1998)). The energy partitioning algorithm is activated either by supplying the ENEPART
directive:
ENEPART
,[epart],[iepart]
or by giving the parameters as options on the command line.
The epart parameter determines the cutoff distance for (intramolecular) bond lengths (in a.u., default 3 a.u.) and is used to automatically determine the individual monomer subunits of the cluster. The iepart parameter enables the energy partitioning, if set to a value larger than zero (default 1). Additionally, if iepart is set to 2, a list of all intermolecular pair energies and their components is printed.
The output section produced by the energy partitioning algorithm will look similar to the following example:
energy partitioning enabled ! centre groups formed for cutoff [au] = 3.00 1 :O1 H11 H12 2 :O2 H21 H22 energy partitioning relative to centre groups: intramolecular correlation: .43752663 exchange dispersion : .00000037 dispersion energy : .00022425 ionic contributions : .00007637
The centre groups correspond to the individual monomers determined for epart=3. In the present example, two water monomers were found. The correlation energy is partitioned into the four components shown above. The exchange dispersion, dispersion and ionic components reflect directly the related intermolecular components of the complex, while the intramolecular correlation contribution to the interaction energy has to be determined by a supermolecular calculation, i.e. by subtracting the (two) corresponding monomer correlation energies from the intramolecular correlation component of the complex given in the output.
Alternatively, the following form can be used:
ENEPART,RMAX=[
r1,r2,r3,…]
and the program will then print the energy contributions of all pairs in the ranges between the given distances (in bohr, enclosed by square brackets, e.g., enepart,rmax=[0,3,5,7,9,11]
). A second list in which the contributions are given as a function of the number of bonds between the pair domains will also be printed.
Split Coulomb operator treatment
Doing it right
The local correlation methods in Molpro employ localized molecular orbitals (LMOs). PipekMezey localization is recommended, but Boys localization is also possible. The virtual orbital space is spanned by nonorthogonal projected atomic orbitals (PAOs). The local character of this basis makes it possible to introduce two distinct approximations: first, excitations are restricted to domains, which are subspaces of (PAOs) that are spatially close to the orbitals from which the electrons are being excited. Secondly, the orbital pairs are classified according to their importance (based on distance or connectivity criteria), and only strong pairs are treated at the highest level (e.g. CCSD). The remaining weak and distant pairs are treated at the LMP2 level, and very distant pairs are neglected. These approximations lead to linear scaling of the computational resources as a function of the molecular size.
Naturally, such approximation can introduce some errors, and therefore the user has to be more careful than with standard black box methods. On the other hand, the loworder scaling makes it possible to treat much larger systems at high levels of theory than it was possible so far.
This section summarizes some important points to remember when performing local correlation calculations.
Basis sets
For numerical reasons, it is useful to eliminate projected core orbitals, since these may have a very small norm. By default, projected core orbitals are eliminated if their norm is smaller than 0.1 (this behaviour can be changed using the DELCOR
and THRCOR
options). For local calculations we recommend the use of generally contracted basis sets, e.g., the correlation consistent ccpVnZ sets of Dunning and coworkers. For these basis sets the core basis functions are uniquely defined, and will always be eliminated if the defaults for DELCOR
and THRCOR
are used.
The correlation consistent basis sets are also recommended for all density fitting calculations, since optimized fitting basis sets are available for each basis.
Symmetry and Orientation
1. Turn off symmetry! Otherwise, you won’t get appropriately localized orbitals (local orbitals will tend to be symmetry equivalent instead of symmetry adapted). Symmetry can be used only if all atoms are symmetry unique. This allows the local treatment of planar molecules in $C_s$ symmetry. But note that neither the multipole program nor the density fitting programs support symmetry at all, so choose always $C_1$ symmetry for DFcalculations or with the MULTP
option.
To turn off symmetry, add NOSYM
or SYMMETRY,NOSYM
before the geometry block, e.g.
nosym geometry={ O1 H1,O1,roh H2,O1,roh,h1,hoh }
2. Use NOORIENT! We recommend to use the NOORIENT
option in the geometry input, to avoid unintended rotations of the molecule when the geometry changes. This is particularly important for geometry optimizations and for domain restarts in calculations of interaction energies (see section intermolecular interactions).
Localization
By default, PipekMezey localization is used and performed automatically in the beginning of a local correlation calculation. Thus
dfhf !HartreeFock with density fitting dflmp2 !LMP2 using the PipekMezey LMOs
is equivalent to
dfhf !HartreeFock with density fitting locali,pipek !Orbital localization using the PipekMezey criterion dflmp2 !LMP2 using the PipekMezey LMOs
Boys localization can be used as well, but in this case the localization must be done beforehand, e.g.
dfhf !HartreeFock with density fitting locali,boys !Orbital localization using the Boys criterion dflmp2 !LMP2 using the Boys LMOs
Poor localization is sometimes an intrinsic problem, in particular for strongly conjugated systems or when diffuse basis sets are used. This is caused by localization tails due to the overlapping diffuse functions. The problem is particularly frequent in calculations of systems with short bonds, e.g., aromatic molecules. It can be avoided using directive
PIPEK,DELETE=
$n$
with $n=1$ or $2$. This means that the contributions of the $n$ most diffuse basis functions of each angular momentum type are ignored in the localization. This often yields much better localized orbitals when diffuse basis sets are used. For augccpVTZ, $n=2$ has been found to work very well, while for augccpVDZ $n$=1
In rare cases it might also happen that the localization procedure does not converge. It is them possible to choose a secondorder NewtonRaphson localization scheme, using the directive
PIPEK,METHOD=2
,[DELETE=
$n$]
This happens for example in LMP2 optimizations of benzene or other highly symmetric aromatics, where the PipekMezey localization has redundant orbital rotations. Instead of PIPEK,METHOD=2
one can also use option CPLSOLVE=2
.
Normally (default) one can use
PIPEK,METHOD=3
,[DELETE=
$n$]
which performs standard PipekMezey iterations, which usually converge quickly.
Orbital domains
The orbital domains are determined automatically using the procedure of Boughton and Pulay, J. Comput. Chem. 14, 736 (1993) and J. Chem. Phys. 104, 6286 (1996). For higher accuracy the domains can be extended, and in this way the canonical result can be systematically approached (cf. Ref. [1] and section extended domains). Details are described in section options for selection of domains.
In most cases, the domain selection is uncritical for saturated molecules. Nevertheless, in particular for delocalized systems, it is recommended always to check the orbital domains, which are printed in the beginning of each local calculation. For such checking, the option DOMONLY=1
can be used to stop the calculation after the domain generation. The orbital domains consist of all basis functions for a subset of atoms. These atoms are selected so that the domain spans the corresponding localized orbital with a preset accuracy (alterable with option THRBP
). A typical domain output, here for water, looks like this:
Orbital domains Orb. Atom Charge Crit. 2.1 1 O1 1.17 0.00 3 H2 0.84 1.00 3.1 1 O1 2.02 1.00 4.1 1 O1 1.96 1.00 5.1 1 O1 1.17 0.00 2 H1 0.84 1.00
This tells you that the domains for orbitals 2.1 and 5.1 comprise the basis functions of the oxygen atom and and one hydrogen atom, while the domains for orbitals 3.1 and 4.1 consist of the basis function on oxygen only. The latter ones correspond to the oxygen lone pairs, the former to the two OH bonds, and so this is exactly what one would expect. For each domain of AOs, corresponding projected atomic orbitals (PAOs) are generated, which span subspaces of the virtual space and into which excitations are made. Options which affect the domain selection are described in section options for selection of domains. Improper domains can result from poorly localized orbitals (see section localization or a forgotten NOSYM
directive. This does not only negatively affect performance and memory requirements, but can also lead to unexpected results.
The default for the selection criterion THRBP
is 0.98. This works usually well for small basis sets like ccpVDZ. For larger basis sets like ccpVTZ we recommend to use a slightly larger value of 0.985 to ensure that enough atoms are included in each domain. For ccpVQZ recommend THRBP=0.990
is recommended. In cases of doubt, compare the domains you get with a smaller basis (e.g., ccpVDZ).
The choice of domains usually has only a weak effect on nearequilibrium properties like equilibrium geometries and harmonic vibrational frequencies. More critical are energy differences like reaction energies or barrier heights. In cases where the electronic structure strongly changes, e.g., when the number of double bonds changes, it is recommended to compare DFLMP2 and DFMP2 results before performing expensive LCCSD(T) calculations. More balanced results and smooth potentials can be obtained using the MERGEDOM
directive, see section domain Merging (MERGEDOM).
Freezing domains
In order to obtain smooth potential energy surfaces, domains must be frozen. The domain information can be stored using the SAVE
option and recovered using the START
option. Alternatively, the SAVE
and START
can be used, see section saving the wavefunction (SAVE). In the latter case, also the CCSD amplitudes are saved/restarted. Freezing domains is particularly important in calculations of intermolecular interactions, see section intermolecular interactions. Domains that are appropriate for larger ranges of geometries, such as reaction pathways, can be generated using the MERGEDOM
directive, section domain Merging (MERGEDOM). The domains are automatically frozen in geometry optimizations and frequency calculations, see section gradients and frequency calculations. Note, however, that this automatic procedure only works if a single local calculation is involved in the optimization. In case of optimizations with counterpoise correction the domains for the complex and each monomer must be frozen individually in different records using the SAVE
and START
directives.
Pair Classes
The strong, close, weak and distant pairs are selected using distance or connectivity criteria as described in more detail in section options for selection of pair classes. Strong pairs are treated by CCSD, all other pairs by LMP2. In triples calculations, all orbital triples $(ijk)$ are included for which $(ij)$, $(ik)$, and $(jk)$ are close pairs. In addition, one of these pairs is restricted to be strong. The triples energy depends on the strong and close pair amplitudes. The close pair amplitudes are taken from the LMP2 calculation. Thus, increasing the distance or connectivity criteria for close and weak pairs will lead to more accurate triples energies. While for near equilibrium properties like geometries and harmonic vibrational frequencies the default values are normally appropriate, larger distance criteria are sometimes needed when computing energy differences, in particular barrier heights. In cases of doubt, RWEAK
should first be increased until convergence is reached, and then RCLOSE
can be varied as well. Such tests can be performed with small basis sets like ccpVDZ, and the optimized values then be used in the final calculations with large basis sets.
Gradients and frequency calculations
Geometry optimizations [1517] and numerical frequency calculations [1820] can be performed using analytical energy gradients [1517] for local MP2. LMP2 geometry optimizations are particularly attractive for weakly bound systems, since virtually BSSE free structures are obtained (see section intermolecular interactions and Refs. [2123]). For reasons of efficiency it is strongly advisable to use the DFLMP2 Gradient [17] for all geometry optimizations. Setting SCSGRD=1
on the DFLMP2
command or DFIT
directive activates the gradient with respect to Grimmes SCS scaled MP2 energy functional (see also section DFIT). Analytical energy gradients are not yet available for the multipole approximation of distant pairs, and therefore MULTP
cannot be used in geometry optimizations or frequency calculations.
In geometry optimizations, the domains are allowed to vary in the initial optimization steps. When the stepsize drops below a certain threshold (default 0.01) the domains are automatically frozen. In numerical Hessian or frequency calculations the domains are also frozen. It is therefore not necessary to include SAVE
and START
options.
Particular care must be taken in optimizations of highly symmetric aromatic systems, like, e.g., benzene. In $D_{6h}$ symmetry, the localization of the $\pi$orbitals is not unique, i.e., the localized orbitals can be rotated around the $C_6$ axis without changing the localization criterion. This redundancy is lost if the symmetry is slightly distorted, which can lead to sudden changes of the localized orbitals. If now the domains are kept fixed using the SAVE
and START
options, a large error in the energy might result. On the other hand, if the domains are not kept fixed, their size and quality might change during the optimization, again leading to spurious energy changes and divergence of the optimization.
The best way to avoid this problem is to use the MERGEDOM=1
option (see section options for selection of domains). If this option is given, the domains for the $\pi$ orbitals will comprise the basis functions of all six carbon atoms, and the energy will be invariant with respect to unitary transformations among the three $\pi$ orbitals. Note that this problem does not occur if the symmetry of the aromatic system is lowered by a substituent.
Redundant orbital rotations can also lead to convergence difficulties of the PipekMezey
localization. This can be overcome by using
PIPEK,METHOD=2
With METHOD=2
, the second derivatives of the localization criterion with respect to the orbital rotations is computed and diagonalized, and rotations corresponding to zero eigenvalues are eliminated. This method converges quadratically and should be used when highly symmetric aromatic systems are optimized. With METHOD=3
(default) standard PipekMezey method are performed.
Finally, we note that the LMP2 gradients are quite sensitive to the accuracy of the SCF convergence (as is also the case for MP2). If very accurate structures are required, or if numerical frequencies are computed from the gradients, the default SCF accuracy might be insufficient. We recommend in such cases to add an ACCU,14
directive (possibly even ACCU,16
) after the HF
command. Indicative of insufficient SCF accuracy are small positive energy changes near the end of the geometry optimization.
Intermolecular interactions
Local methods are particularly useful for the calculation of weak intermolecular interactions since the basis set superposition error (BSSE) is largely reduced [1,13,14] and counterpoise corrections are usually not necessary (provided the BSSE of the underlying HartreeFock is small). However, one must be careful to define the domains properly and to include all intermolecular pairs at the highest computational level. A convenient way to define appropriate domains and pair lists is to use the option INTERACT=1
. If this option is given, individual molecules are identified automatically and all intermolecular pairs are automatically treated as strong pairs and included in the LCCSD. Similarly, appropriate triples lists are generated for LCCSD(T) calculations. It is required that all orbital domains are located on individual molecules. Note however that the inclusion of the intermolecular pairs strongly increases the number of strong pairs and triples, and therefore highlevel calculations can become very expensive.
For calculations of interaction potentials of weakly interacting systems, the domains of the subsystems should be determined at a very large distance and saved using the SAVE=
record option on the LOCAL
or MULTP
directive, or the SAVE
directive (see section saving the wavefunction (SAVE)). If the asymptotic energy is not needed it is sufficient to do this initial calculation using option DOMONLY=1
). These domains should then be reused in the subsequent calculations at all other intermolecular distances by using the START
=record option or the START
directive (see section restarting a calculation (START)). Only in this way the basis set superposition error is minimized and normally negligible (of course, this does not affect the BSSE for the SCF, and therefore the basis set should be sufficiently large to make the SCF BSSE negligible).
Usually, diffuse basis functions are important for obtaining accurate intermolecular interactions. When diffuse basis sets are used, it may happen that the PipekMezey localization does not yield well localized orbitals. This problem can in most cases be overcome by using the directive
PIPEK,DELETE=
$n$
as described in section localization
A final warning concerns local density fitting (see sections Densityfitted LMP2 (DFLMP2), LDCSD (DFLDCSD), and coupled cluster (DFLCCSD(T0)) and density fitting): local fitting must not be used in counterpoise calculations, since no fitting functions would be present on the dummy atoms and this can lead to large errors.
For examples and discussions of these aspects see Refs. [2123]
Efficient and accurate treatment of close pairs in local CCSD(T)
As discussed above, efficient LCCSD(T) calculations employ a hierarchical tretament of the LMO pairs, which are distinguished into the pair classes strong, close, weak, and very distant. Only the strong pairs are actually treated at the full coupled cluster level, while for close and weak pairs only LMP2 is used (and very distant pairs are omitted altogether). Consequently, intermolecular pairs in bigger molecules or clusters, which are usually not strong, are only treated at the LMP2 level, which is often poor (since van der Waals interactions are only captured at the lowest level of perturbation theory) and compromises the otherwise high accuracy of coupled cluster like LCCSD(T). By setting the keyword how_treatclswk
on the local
card appropriately, it is possible to treat the close pairs at a higher level than LMP2, i.e.,
how_treatclswk=1
local MP2how_treatclswk=2
direct local RPAhow_treatclswk=3
local ringCCD3how_treatclswk=4
local ringCCDhow_treatclswk=5
local CCD[S]R$^{6}$
These methods for close pairs are discussed in detail in Refs. [2425]. The philosophy is to include all diagrams with R$^{6}$ decay behavior up to a certain order of perturbation theory. E.g., LringCCD3 denotes local ring CCD up to thirdorder diagrams, and LCCD[S]R$^{6}$ is full local CCD with a perturbative singles correction, yet omitting all diagrams with decay behavior quicker than R$^{6}$ (like ladder diagrams). As shown in Ref. [25], the LCCD[S]R$^{6}$ provides excellent intermolecular interaction energies, virtually indistinguishable from LCCSD, at a much lower cost (mainly due to the omission of ladders). We therefore recommend to use either LringCCD3 or LCCD[S]R$^{6}$. Note that it is essential to couple strong and close pairs in such calculations, i.e., to set keepcls=1
(vide supra). As an example, the input segment to set all intermolecular pairs to close and treat them at the level of LCCD[S]R$^{6}$ is
local,keepcls=1,interact=1,interpair=1,how_treatclswk=5
See also the input example h2o_dimer_lccsdhybrids.inp.
Densityfitted LMP2 (DFLMP2), LDCSD (DFLDCSD), and coupled cluster (DFLCCSD(T0))
[sec:locdf]
Densityfitting LMP2, LDCSD and LCCSD calculations can be performed by adding the prefix DF
to the command name. The input is as follows:
DFLMP2
,[options]
DFLCCSD(T)
,[options]
Options for density fitting can be mixed with any options for LOCAL
. Options for density fitting can also be given on a DFIT
directive (see section density fitting).
The most important options for density fitting in local methods are
BASIS_MP2
=string Fitting basis set used in LMP2 and in LCCSD for integrals with up to 2 external orbitals. If a correlation consistent basis set is used (e.g. ccpVTZ) the corresponding fitting basis for MP2 us used by default (ccpVTZ/MP2FIT). Otherwise the fitting basis set must be defined in a preceding basis block (see section basis input).BASIS_CCSD
=string Fitting basis set used in LCCSD for integrals over 3 and 4external orbitals. The default isBASIS_MP2
and this is usually sufficient. However, the accurate approximation of 4external integrals in LCCSD requires larger fitting basis sets than LMP2. Therefore, in order to minimize fitting errors, it is recommended to use the next larger fitting basis, e.g.,BASIS_CCSD=VQZ
for orbital basisVTZ
.LOCFIT
=value: IfLOCFIT=1
local fitting is enabled. This is necessary to achieve linear scaling in DFLMP2 (see Refs. [1114]). The errors introduced by local fitting are usually very small, but there are some exceptions. For instance,LOCFIT=1
must not be used in counterpoise calculations, see section intermolecular interactions) Note that for small moleculesLOCFIT=1
can be more expensive thanLOCFIT=0
.
For further details and options for density fitting see section density fitting.