1% 2% $Id$ 3% 4\label{sec:dft} 5 6The NWChem density functional theory (DFT) module uses the 7Gaussian basis set approach to compute 8closed shell and open shell densities and Kohn-Sham orbitals 9in the: 10\begin{itemize} 11\item local density approximation (LDA), 12\item non-local density approximation (NLDA), 13\item local spin-density approximation (LSD), 14\item non-local spin-density approximation (NLSD), 15\item non-local meta-GGA approximation (metaGGA), 16\item any empirical mixture of local and non-local approximations 17(including exact exchange), and 18\item asymptotically corrected exchange-correlation potentials. 19\end{itemize} 20 21The formal scaling of the DFT computation can be reduced by choosing 22to use auxiliary Gaussian basis sets to fit the charge density (CD) and/or 23fit the exchange-correlation (XC) potential. 24 25DFT input is provided using the compound \verb+DFT+ directive 26\begin{verbatim} 27 DFT 28 ... 29 END 30\end{verbatim} 31The actual DFT calculation will be performed when the input module 32encounters the \verb+TASK+ directive (Section \ref{sec:task}). 33\begin{verbatim} 34 TASK DFT 35\end{verbatim} 36 37Once a user has specified a geometry and a Kohn-Sham orbital basis set 38the DFT module can be invoked with no input directives (defaults 39invoked throughout). There are sub-directives which allow for 40customized application; those currently provided as options for 41the DFT module are: 42\begin{verbatim} 43 VECTORS [[input] (<string input_movecs default atomic>) || \ 44 (project <string basisname> <string filename>)] \ 45 [swap [alpha||beta] <integer vec1 vec2> ...] \ 46 [output <string output_filename default input_movecs>] \ 47 48 49 XC [[acm] [b3lyp] [beckehandh] [pbe0]\ 50 [becke97] [becke97-1] [becke97-2] [becke97-3] [becke97-d] [becke98] \ 51 [hcth] [hcth120] [hcth147]\ 52 [hcth407] [becke97gga1] [hcth407p]\ 53 [mpw91] [mpw1k] [xft97] [cft97] [ft97] [op] [bop] [pbeop]\ 54 [xpkzb99] [cpkzb99] [xtpss03] [ctpss03] [xctpssh]\ 55 [b1b95] [bb1k] [mpw1b95] [mpwb1k] [pw6b95] [pwb6k] [m05] [m05-2x] [vs98] \ 56 [m06] [m06-hf] [m06-L] [m06-2x] \ 57 [HFexch <real prefactor default 1.0>] \ 58 [becke88 [nonlocal] <real prefactor default 1.0>] \ 59 [xperdew91 [nonlocal] <real prefactor default 1.0>] \ 60 [xpbe96 [nonlocal] <real prefactor default 1.0>] \ 61 [gill96 [nonlocal] <real prefactor default 1.0>] \ 62 [lyp <real prefactor default 1.0>] \ 63 [perdew81 <real prefactor default 1.0>] \ 64 [perdew86 [nonlocal] <real prefactor default 1.0>] \ 65 [perdew91 [nonlocal] <real prefactor default 1.0>] \ 66 [cpbe96 [nonlocal] <real prefactor default 1.0>] \ 67 [pw91lda <real prefactor default 1.0>] \ 68 [slater <real prefactor default 1.0>] \ 69 [vwn_1 <real prefactor default 1.0>] \ 70 [vwn_2 <real prefactor default 1.0>] \ 71 [vwn_3 <real prefactor default 1.0>] \ 72 [vwn_4 <real prefactor default 1.0>] \ 73 [vwn_5 <real prefactor default 1.0>] \ 74 [vwn_1_rpa <real prefactor default 1.0>] \ 75 [xtpss03 [nonlocal] <real prefactor default 1.0>] \ 76 [ctpss03 [nonlocal] <real prefactor default 1.0>] \ 77 [bc95 [nonlocal] <real prefactor default 1.0>] \ 78 [xpw6b95 [nonlocal] <real prefactor default 1.0>] \ 79 [xpwb6k [nonlocal] <real prefactor default 1.0>] \ 80 [xm05 [nonlocal] <real prefactor default 1.0>] \ 81 [xm05-2x [nonlocal] <real prefactor default 1.0>] \ 82 [cpw6b95 [nonlocal] <real prefactor default 1.0>] \ 83 [cpwb6k [nonlocal] <real prefactor default 1.0>] \ 84 [cm05 [nonlocal] <real prefactor default 1.0>] \ 85 [cm05-2x [nonlocal] <real prefactor default 1.0>]] \ 86 [xvs98 [nonlocal] <real prefactor default 1.0>]] \ 87 [cvs98 [nonlocal] <real prefactor default 1.0>]] \ 88 [xm06-L [nonlocal] <real prefactor default 1.0>]] \ 89 [xm06-hf [nonlocal] <real prefactor default 1.0>]] \ 90 [xm06 [nonlocal] <real prefactor default 1.0>]] \ 91 [xm06-2x [nonlocal] <real prefactor default 1.0>]] \ 92 [cm06-L [nonlocal] <real prefactor default 1.0>]] \ 93 [cm06-hf [nonlocal] <real prefactor default 1.0>]] \ 94 [cm06 [nonlocal] <real prefactor default 1.0>]] \ 95 [cm06-2x [nonlocal] <real prefactor default 1.0>]] 96 97 98 99 CONVERGENCE [[energy <real energy default 1e-7>] \ 100 [density <real density default 1e-5>] \ 101 [gradient <real gradient default 5e-4>] \ 102 [dampon <real dampon default 0.0>] \ 103 [dampoff <real dampoff default 0.0>] \ 104 [diison <real diison default 0.0>] \ 105 [diisoff <real diisoff default 0.0>] \ 106 [levlon <real levlon default 0.0>] \ 107 [levloff <real levloff default 0.0>] \ 108 [ncydp <integer ncydp default 2>] \ 109 [ncyds <integer ncyds default 30>] \ 110 [ncysh <integer ncysh default 30>] \ 111 [damp <integer ndamp default 0>] [nodamping] \ 112 [diis [nfock <integer nfock default 10>]] \ 113 [nodiis] [lshift <real lshift default 0.5>] \ 114 [nolevelshifting] \ 115 [hl_tol <real hl_tol default 0.1>] \ 116 [rabuck [n_rabuck <integer n_rabuck default 25>]] 117 118 119 GRID [(xcoarse||coarse||medium||fine||xfine) default medium] \ 120 [(gausleg||lebedev ) default lebedev ] \ 121 [(becke||erf1||erf2||ssf) default erf1] \ 122 [(euler||mura||treutler) default mura] \ 123 [rm <real rm default 2.0>] \ 124 [nodisk] 125 126 127 TOLERANCES [[tight] [tol_rho <real tol_rho default 1e-10>] \ 128 [accCoul <integer accCoul default 8>] \ 129 [radius <real radius default 25.0>]] 130 131 132 [(LB94||CS00 <real shift default none>)] 133 134 DECOMP 135 ODFT 136 DIRECT 137 SEMIDIRECT [filesize <integer filesize default disksize>] 138 [memsize <integer memsize default available>] 139 [filename <string filename default $file_prefix.aoints$>] 140 INCORE 141 ITERATIONS <integer iterations default 30> 142 MAX_OVL 143 MULLIKEN 144 DISP 145 MULT <integer mult default 1> 146 NOIO 147 PRINT||NOPRINT 148\end{verbatim} 149% [accCDfunc <integer accAOfunc default 20>] \ 150% [store_wght] [nquad_task <integer nquad_task default 1>] \ 151 152The following 153sections describe these keywords and 154optional sub-directives that can be specified for a \verb+DFT+ calculation 155in NWChem. 156 157\section{Specification of Basis Sets for the DFT Module} 158 159The DFT module requires at a minimum the basis set for the Kohn-Sham 160molecular orbitals. This basis set must be in the default basis set named 161{\tt "ao basis"}, or it must be assigned to this default name using the 162\verb+SET+ directive (see Section \ref{sec:set}). 163 164In addition to the basis set for the Kohn-Sham orbitals, 165the charge density fitting basis set can also be specified in the 166input directives for the DFT module. This basis set is used for the 167evaluation of the Coulomb potential in the Dunlap scheme\footnote{B.I.~Dunlap, 168J.W.D.~Connolly and J.R.~Sabin, J.~Chem.~Phys.~{\bf 71}, 4993 (1979)}. 169The charge density fitting basis set must have the name {\tt "cd basis"}. 170This can be the actual name of a basis set, or a basis set can be 171assigned this name using the \verb+SET+ directive, as described in 172Section \ref{sec:set}. If this basis set is not defined by input, 173the $O(N^4)$ exact Coulomb contribution is computed. 174 175The user also has the option of specifying a third basis set for the 176evaluation of the exchange-correlation potential. This basis set must 177have the name {\tt "xc basis"}. If this basis set is not specified 178by input, the exchange contribution (XC) is evaluated by numerical 179quadrature. In most applications, this approach is efficient enough, 180so the {\tt "xc basis"} basis set is not generally required. 181 182For the DFT module, the input options for defining the basis sets in a given 183calculation can be summarized as follows; 184\begin{itemize} 185\item {\tt "ao basis"} -- Kohn-Sham molecular orbitals; required for all 186calculations 187\item {\tt "cd basis"} -- charge density fitting basis set; optional, but 188recommended for evaluation of the Coulomb potential 189\item {\tt "xc basis"} -- exchange-correlation (XC) fitting basis set; 190optional, and usually not recommended 191\end{itemize} 192 193 194\section{{\tt VECTORS} and {\tt MAX\_OVL} --- KS-MO Vectors} 195 196The \verb+VECTORS+ directive is the same as that in the SCF module 197(Section \ref{sec:vectors}). Currently, the \verb+LOCK+ keyword 198is not supported by the DFT module, however the directive 199\begin{verbatim} 200 MAX_OVL 201\end{verbatim} 202has the same effect. 203 204\section{{\tt XC} and {\tt DECOMP} --- Exchange-Correlation Potentials} 205\label{sec:xc} 206\begin{verbatim} 207 XC [[acm] [b3lyp] [beckehandh] [pbe0]\ 208 [becke97] [becke97-1] [becke97-2] [becke97-3] [becke98] [hcth] [hcth120] [hcth147] \ 209 [hcth407] [becke97gga1] [hcth407p] \ 210 [optx] [hcthp14] [mpw91] [mpw1k] [xft97] [cft97] [ft97] [op] [bop] [pbeop]\ 211 [HFexch <real prefactor default 1.0>] \ 212 [becke88 [nonlocal] <real prefactor default 1.0>] \ 213 [xperdew91 [nonlocal] <real prefactor default 1.0>] \ 214 [xpbe96 [nonlocal] <real prefactor default 1.0>] \ 215 [gill96 [nonlocal] <real prefactor default 1.0>] \ 216 [lyp <real prefactor default 1.0>] \ 217 [perdew81 <real prefactor default 1.0>] \ 218 [perdew86 [nonlocal] <real prefactor default 1.0>] \ 219 [perdew91 [nonlocal] <real prefactor default 1.0>] \ 220 [cpbe96 [nonlocal] <real prefactor default 1.0>] \ 221 [pw91lda <real prefactor default 1.0>] \ 222 [slater <real prefactor default 1.0>] \ 223 [vwn_1 <real prefactor default 1.0>] \ 224 [vwn_2 <real prefactor default 1.0>] \ 225 [vwn_3 <real prefactor default 1.0>] \ 226 [vwn_4 <real prefactor default 1.0>] \ 227 [vwn_5 <real prefactor default 1.0>] \ 228 [vwn_1_rpa <real prefactor default 1.0>]] 229\end{verbatim} 230 231The user has the option of specifying the exchange-correlation 232treatment in the DFT Module (see table \ref{tablexc}). 233 The default exchange-correlation 234functional is defined as the local density approximation (LDA) for 235closed shell systems and its counterpart the local spin-density (LSD) 236approximation for open shell systems. Within this approximation the 237exchange functional is the Slater $\rho^{1/3}$ functional (from 238J.C.~Slater, {\sl Quantum Theory of Molecules and Solids, Vol.~4: The 239 Self-Consistent Field for Molecules and Solids} (McGraw-Hill, New 240York, 1974)), and the correlation functional is the Vosko-Wilk-Nusair 241(VWN) functional (functional V) (S.J.~Vosko, L.~Wilk and M.~Nusair, 242Can.~J.~Phys.~{\bf 58}, 1200 (1980)). The parameters used in this 243formula are obtained by fitting to the Ceperley and 244Alder\footnote{D.M.~Ceperley and B.J.~Alder, Phys. Rev. Lett. {\bf 245 45}, 566 (1980).} 246Quantum Monte-Carlo solution of the {\em 247 homogeneous electron gas}. 248 249These defaults can be invoked explicitly by specifying the following 250keywords within the DFT module input directive, \verb+XC slater vwn_5+. 251 252 253That is, this statement in the input file 254\begin{verbatim} 255dft 256 XC slater vwn_5 257end 258task dft 259\end{verbatim} 260 261is equivalent to the simple line 262\begin{verbatim} 263task dft 264\end{verbatim} 265 266 267The \verb+DECOMP+ directive causes the components of the energy 268corresponding to each functional to be printed, rather than just the 269total exchange-correlation energy which is the default. You can see 270an example of this directive in the sample input in 271Section \ref{sec:DFTsample}. 272 273 274Many alternative exchange and correlation functionals are available to 275the user as listed in table \ref{tablexc}. The following sections describe 276how to use these options. 277 278\subsection{Exchange-Correlation Functionals} 279 280There are several Exchange and Correlation functionals in addition to the 281default {\tt slater} and {\tt vwn\_5} 282functionals. These are either local or gradient-corrected functionals (GCA); 283a full list can be found in table \ref{tablexc}. 284 285%\begin{itemize} 286% 287%\item the Becke 1998 gradient-corrected functional (see A.D.~Becke, 288% Phys. Rev. A {\bf 38}, 3098 (1988)), is invoked by specifying 289% 290%\begin{verbatim} 291% XC becke88 292%\end{verbatim} 293% 294%\item the Perdew 1991 gradient-corrected exchange functional (J.P. Perdew, 295% J.A.~Chevary, S.H.~Vosko, K.A.~Jackson, M.R.~Pederson, D.J.~Singh 296% and C.~Fiolhais, Phys. Rev. B {\bf 46}, 6671 (1992)), is invoked by specifying 297% 298%\begin{verbatim} 299% XC xperdew91 300%\end{verbatim} 301% 302%\item the Perdew-Burke-Ernzerhof gradient-corrected exchange functional \\ 303% (J.P.~Perdew, K.~Burke and M.~Ernzerhof, Physical Review Letters 304%{\bf 77}, 3865 (1996); {\bf 78}, 1396 (1997))), is invoked by specifying 305% 306%\begin{verbatim} 307% XC xpbe96 308%\end{verbatim} 309% 310%\item the Gill 1996 gradient-corrected exchange functional (P.W.Gill , Mol. Phys. {\bf 89}, 433 (1996)), is invoked by specifying 311% 312%\begin{verbatim} 313% XC gill96 314%\end{verbatim} 315 316 317%\item 318The Hartree-Fock exact exchange functional, (which has $O(N^4)$ 319computation expense), is invoked by specifying 320\begin{verbatim} 321 XC HFexch 322\end{verbatim} 323%\end{itemize} 324 325Note that the user also has the ability to include only the local or 326nonlocal contributions of a given functional. In addition the user 327can specify a multiplicative prefactor (the variable 328\verb+<prefactor>+ in the input) for the local/nonlocal component or 329total. An example of this might be, 330\begin{verbatim} 331 XC becke88 nonlocal 0.72 332\end{verbatim} 333The user should be aware that the Becke88 local component is simply 334the Slater exchange and should be input as such. 335 336Any combination of the supported exchange functional options can be 337used. For example the popular Gaussian B3 exchange could be specified 338as: 339\begin{verbatim} 340 XC slater 0.8 becke88 nonlocal 0.72 HFexch 0.2 341\end{verbatim} 342 343 344%\subsection{Correlation Functionals} 345 346%In addition to the default \verb+vwn_5+ correlation functional, the user has 347%alternative correlation functionals as listed in table \ref{tablexc}. 348%to choose from: lyp, perdew81, 349%perdew86, perdew91, pw91lda, \verb+vwn_1+, \verb+vwn_2+, \verb+vwn_3+, 350%\verb+vwn_4+, and \verb+vwn_1_rpa+. 351 352%As in the exchange functional input, individual local/nonlocal 353%components as well as multiplicative prefactors can be invoked where 354%appropriate. Each of the correlation functionals is listed below along with 355%appropriate citation in table \ref{tablexc}. 356 357\sloppy 358 359%\begin{itemize} 360%\item VWN local density functionals; S.J.~Vosko, L.~Wilk and M.~Nusair, 361% Can.~J.~Phys.~{\bf 58}, 1200 (1980); all five (5) functionals as 362% described in this paper (addressed in the paper as I - V) have been 363% implemented. These functionals can be invoked with the keywords: 364%\begin{verbatim} 365% XC vwn_1 366% XC vwn_2 367% XC vwn_3 368% XC vwn_4 369% XC vwn_5 370%\end{verbatim} 371% 372% Note that functionals; \verb+vwn_2+, \verb+vwn_3+, and \verb+vwn_4+ 373% require both sets of parameters (the Monte Carlo parameters of 374% Ceperley and Alder and VWN's RPA parameters) used in fitting the 375% homogeneous electron gas correlation energy. Functionals 376% \verb+vwn_1+ and \verb+vwn_5+ require only the Monte Carlo fitting 377% parameters. In order to reproduce results in the literature another 378% functional was added; the \verb+vwn_1_rpa+. This is the original 379% \verb+vwn_1+ functional with RPA parameters as opposed to the 380% prescribed Monte Carlo parameters. This functional can be invoked 381% with the keyword, 382%\begin{verbatim} 383% XC vwn_1_rpa 384%\end{verbatim} 385% 386%\item Perdew81 local density functional; J.~P.~Perdew and A.~Zunger, 387% Phys.~Rev.~B {\bf23}, 5048 (1981). This functional can be invoked with the 388% keyword, 389%\begin{verbatim} 390% XC perdew81 391%\end{verbatim} 392% 393%\item Perdew \& Wang 1991 local density functional; J.P.~Perdew 394% and Y.~Wang, Phys. Rev. B {\bf 45}, 13244 (1992). The parameters 395% used in this formula are obtained by fitting to the Ceperley and 396% Alder Quantum Monte Carlo solution of the {\em 397% homogeneous electron gas}. This functional can be invoked with the 398% keyword, 399%\begin{verbatim} 400% XC pw91lda 401%\end{verbatim} 402% 403%\item Perdew86 gradient-corrected functional; J.~P.~Perdew, Phys.~Rev.~B 404% {\bf33}, 8822 (1986). Note that this is a nonlocal functional and 405% in the absence of any local functional specification the local 406% component is defaulted to the perdew81 local correlation 407% functional. This functional can be invoked with the 408% keyword, 409%\begin{verbatim} 410% XC perdew86 411%\end{verbatim} 412% 413%\item Perdew91 gradient-corrected functional; J.P.~Perdew, 414% J.A.~Chevary, S.H.~Vosko, K.A.~Jackson, M.R.~Pederson, D.J.~Singh 415% and C.~Fiolhais, Phys. Rev. B {\bf 46}, 6671 (1992). Note that this 416% is a nonlocal functional and in the absence of any local functional 417% specification the local component is defaulted to the \verb+pw91lda+ local 418% correlation functional. This functional can be invoked with the keyword, 419%\begin{verbatim} 420% XC perdew91 421%\end{verbatim} 422% 423%\item the Perdew-Burke-Ernzerhof gradient-corrected correlation functional \\ 424% (J.P.~Perdew, K.~Burke and M.~Ernzerhof, Physical Review Letters 425%{\bf 77}, 3865 (1996); {\bf 78}, 1396 (1997))), is invoked by specifying 426% 427%\begin{verbatim} 428% XC cpbe96 429%\end{verbatim} 430 431%\item LYP gradient-corrected functional; C.~Lee, W.~Yang and 432% R.~G.~Parr, Phys.~Rev.~B {\bf 37}, 785 (1988). Note that this 433% is a local and nonlocal functional but cannot be conveniently split 434% into the individual components. The option to scale the total remains. 435% This functional can be invoked with the keyword, 436%\begin{verbatim} 437% XC lyp 438%\end{verbatim} 439% 440%\end{itemize} 441 442\fussy 443 444Any combination of the supported correlation functional options can be 445used. For example B3LYP could be specified as: 446\begin{verbatim} 447XC vwn_1_rpa 0.19 lyp 0.81 HFexch 0.20 slater 0.80 becke88 nonlocal 0.72 448\end{verbatim} 449and X3LYP as: 450\begin{verbatim} 451xc vwn_1_rpa 0.129 lyp 0.871 hfexch 0.218 slater 0.782 \ 452becke88 nonlocal 0.542 xperdew91 nonlocal 0.167 453\end{verbatim} 454 455 456\subsection{Combined Exchange and Correlation Functionals} 457 458In addition to the options listed above for the exchange and correlation 459functionals, the user has the alternative of specifying combined exchange and 460correlation functionals. A complete list of the available functionals 461appears in table \ref{tablexc}. 462 463 The available hybrid functionals 464(where a Hartree-Fock Exchange component is present) consist of the Becke 465``{\it half and half}'' (see A.D.~Becke, J.~Chem.~Phys.~98, 1372 (1992)), the 466adiabatic connection method (see A.D.~Becke, J.~Chem.~Phys.~98, 5648 467(1993)), B3LYP (popularized by Gaussian9X), Becke 1997 468(``Becke V'' paper: A.D.Becke, J. Chem. Phys., {\bf 107}, 8554 (1997)). 469 470%These options can be invoked by specifying any of the following input lines, 471%\begin{verbatim} 472% XC beckehandh 473% XC acm 474% XC b3lyp 475% XC becke97 476% Xc becke97-1 477% xc becke98 478% xc pbe0 479% xc hcth 480% xc hcth120 481% xc hcth147 482% xc hcth407 483% xc hcthp14 484% xc becke97gga1 485%\end{verbatim} 486 487The keyword \verb+beckehandh+ specifies that the exchange-correlation energy will be 488computed as 489\begin{eqnarray*} 490E_{XC} \ \approx \ \frac{1}{2} E^{\rm HF}_X + \frac{1}{2} E^{\rm 491 Slater}_{X} + \frac{1}{2} E^{\rm PW91LDA}_{C} 492\end{eqnarray*} 493We know this is NOT the correct Becke prescribed implementation which 494requires the XC potential in the energy expression. But this is what 495is currently implemented as an approximation to it. 496 497%\clearpage 498% 499 500The keyword \verb+acm+ specifies that the exchange-correlation energy 501is computed as 502\begin{eqnarray*} 503E_{XC} \ &=& \ a_0 E^{\rm HF}_X + (1-a_0) E^{\rm Slater}_{X} + 504a_X \Delta E^{\rm Becke88}_{X} + E^{\rm VWN}_C + a_C \Delta E^{Perdew91}_C \\ 505& &{\rm where } \\ 506a_0 &=& 0.20, \ a_X = 0.72, \ a_C = 0.81 507\end{eqnarray*} 508and $\Delta$ stands for a non-local component. 509 510 511The keyword \verb+b3lyp+ specifies that the exchange-correlation energy 512is computed as 513\begin{eqnarray*} 514E_{XC} \ &=& \ a_0 E^{\rm HF}_X + (1-a_0) E^{\rm Slater}_{X} + 515a_X \Delta E^{\rm Becke88}_{X} + (1-a_C)E^{\rm VWN\_1\_RPA}_C + a_C E^{LYP}_C \\ 516& &{\rm where } \\ 517a_0 &=& 0.20, \ a_X = 0.72, \ a_C = 0.81 518\end{eqnarray*} 519 520 521%The keyword \verb+becke97-1+ specifies the hybrid exchange-correlation energy 522%derived by Handy et al by re-fitting the Becke 1997 functional 523%(F.A.Hamprecht, A.J.Cohen, D.J.Tozer and N.C.Handy, 524% J. Chem. Phys. {\bf 109}, 6264 (1998)) 525% 526%The keyword \verb+hcth+ specifies the exchange-correlation energy 527%functional derived by Hamprecht-Cohen-Tozer-Handy 528%(this is not a hybrid functional) 529%(F.A.Hamprecht, A.J.Cohen, D.J.Tozer and N.C.Handy, 530% J. Chem. Phys. {\bf 109}, 6264 (1998)) 531 532 533 534\subsection{Meta-GGA Functionals} 535 536 537One way to calculate meta-GGA energies is to use 538 orbitals and densities 539from fully self-consistent GGA or LDA calculations 540and run them in one iteration in the meta-GGA functional. 541It is expected that meta-GGA energies obtained 542this way will be close to fully self consistent 543meta-GGA calculations. 544 545It is possible to calculate metaGGA energies both ways 546in NWChem, that is, self-consistently or with 547GGA/LDA orbitals and densities. However, since second derivatives 548are not available for metaGGAs, in order to 549calculate frequencies, one must use 550\verb+ task dft freq numerical+. 551A sample file with this is shown below, 552in \ref{sec:DFTsample}. In this instance, the 553energy is calculated self-consistently and geometry is optimized using the analytical gradients. 554 555 556(For more information on metaGGAs, see 557S. Kurth, J. Perdew, P. Blaha, Int. J. Quant. Chem 75, 889 (1999) 558for a brief description of meta-GGAs, and citations 14-27 559therein for thorough background ) 560 561Note: both TPSS and PKZB correlation 562require the PBE GGA CORRELATION (which is itself dependent on an LDA). 563The decision has been made to 564use these functionals with the accompanying local 565PW91LDA. The user does not have the ability to set 566the local part of these metaGGA functionals. 567 568 569\section{{\tt LB94} and {\tt CS00} --- Asymptotic correction} 570 571The keyword \verb+LB94+ will correct the asymptotic region of 572the \verb+XC+ definition of exchange-correlation {\it potential} by 573the van-Leeuwen--Baerends exchange-correlation {\it potential} that 574has the correct $-1/r$ asymptotic behavior. The total energy will be computed by the 575\verb+XC+ definition of exchange-correlation functional. This scheme is known to 576tend to overcorrect the deficiency of most uncorrected exchange-correlation potentials. 577 578The keyword \verb+CS00+, when supplied with a real value of shift (in atomic units), 579will perform Casida--Salahub '00 asymptotic correction. This is primarily intended 580for use in conjunction with TDDFT and the background of this method is given in more 581detail in Chapter 14. The shift is normally positive (which means that the original 582uncorrected exchange-correlation potential must be shifted down). 583 584When the keyword \verb+CS00+ is specified without the value of shift, the program will 585automatically supply it according to the semi-empirical formula of Zhan, Nichols, and 586Dixon (again, see Chapter 14 for more details and references). As the Zhan's formula 587is calibrated against B3LYP results, it is most meaningful to use this in conjunction 588with the B3LYP functional, although the program does not prohibit (or even warn) the use 589of any other functional. 590 591Sample input files of asymptotically corrected TDDFT calculations can be found in 592Chapter 14. 593 594\section{Sample input file} 595\label{sec:DFTsample} 596A simple example calculates the geometry 597 of water, using the metaGGA functionals 598\verb+ xtpss03+ and \verb+ctpss03+. 599This also highlights some of 600the print features in the DFT module. Note 601that you must use the line 602\verb+task dft freq numerical+ because 603analytic hessians are not available for the metaGGAs: 604\begin{verbatim} 605title "WATER 6-311G* meta-GGA XC geometry" 606echo 607geometry units angstroms 608 O 0.0 0.0 0.0 609 H 0.0 0.0 1.0 610 H 0.0 1.0 0.0 611end 612basis 613 H library 6-311G* 614 O library 6-311G* 615end 616dft 617 iterations 50 618 print kinetic_energy 619 xc xtpss03 ctpss03 620 decomp 621end 622task dft optimize 623 624task dft freq numerical 625\end{verbatim} 626 627 628 629 630%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 631\twocolumn 632\begin{table}[htp] 633\centering 634{\scriptsize 635\begin{tabular}{|l|p{0.12cm}p{0.12cm}cccp{0.35cm}|c|} 636\hline 637Keyword & X & C & GGA & Meta & Hybr. & 2nd & Ref.\\ 638\hline 639{\tt slater }& $\star$ & & & & & Y &[1] \\ 640\hline 641{\tt vwn\_1 }& & $\star$ & & & & Y &[2] \\ 642{\tt vwn\_2 }& & $\star$ & & & & Y &[2] \\ 643{\tt vwn\_3 }& & $\star$ & & & & Y &[2] \\ 644{\tt vwn\_4 }& & $\star$ & & & & Y &[2] \\ 645{\tt vwn\_5 }& & $\star$ & & & & Y &[2] \\ 646{\tt vwn\_1\_rpa}& & $\star$ & & & & Y &[2] \\ 647{\tt perdew81 }& & $\star$ & & & & Y &[3] \\ 648{\tt pw91lda }& & $\star$ & & & & Y &[4] \\ 649\hline 650{\tt becke88 }& $\star$ & & $\star$ & & & Y &[5]\\ 651{\tt xperdew91}& $\star$ & & $\star$ & & & Y &[6]\\ 652{\tt xpbe96 }& $\star$ & & $\star$ & & & Y &[7]\\ 653{\tt gill96 }& $\star$ & & $\star$ & & & Y &[8]\\ 654{\tt optx }& $\star$ & & $\star$ & & & N &[20]\\ 655{\tt mpw91 }& $\star$ & & $\star$ & & & Y &[23]\\ 656{\tt xft97 }& $\star$ & & $\star$ & & & N &[24]\\ 657{\tt rpbe }& $\star$ & & $\star$ & & & Y &[33]\\ 658{\tt revpbe }& $\star$ & & $\star$ & & & Y &[34]\\ 659{\tt xpw6b95 }& $\star$ & & $\star$ & & & N &[36]\\ 660{\tt xpwb6k }& $\star$ & & $\star$ & & & N &[36]\\ 661\hline 662{\tt perdew86 }& & $\star$ & $\star$ & & & Y &[9]\\ 663{\tt lyp }& & $\star$ & $\star$ & & & Y &[10]\\ 664{\tt perdew91 }& & $\star$ & $\star$ & & & Y &[6]\\ 665{\tt cpbe96 }& & $\star$ & $\star$ & & & Y &[7]\\ 666{\tt cft97 }& & $\star$ & $\star$ & & & N &[24]\\ 667{\tt op }& & $\star$ & $\star$ & & & N &[31]\\ 668\hline 669{\tt hcth }& $\star$ & $\star$ & $\star$ & & & N &[11]\\ 670{\tt hcth120 }& $\star$ & $\star$ & $\star$ & & & N &[12]\\ 671{\tt hcth147 }& $\star$ & $\star$ & $\star$ & & & N &[12]\\ 672{\tt hcth407 }& $\star$ & $\star$ & $\star$ & & & N &[19]\\ 673{\tt becke97gga1}& $\star$ & $\star$ & $\star$ & & & N &[18]\\ 674{\tt hcthp14 }& $\star$ & $\star$ & $\star$ & & & N &[21]\\ 675{\tt ft97 }& $\star$ & $\star$ & $\star$ & & & N &[24]\\ 676{\tt htch407p }& $\star$ & $\star$ & $\star$ & & & N &[27]\\ 677{\tt bop }& $\star$ & $\star$ & $\star$ & & & N &[31]\\ 678{\tt pbeop }& $\star$ & $\star$ & $\star$ & & & N &[32]\\ 679\hline 680{\tt xpkzb99 }& $\star$ & & &$\star$& & N &[26]\\ 681{\tt cpkzb99 }& & $\star$ & &$\star$& & N &[26]\\ 682{\tt xtpss03 }& $\star$ & & &$\star$& & N &[28]\\ 683{\tt ctpss03 }& & $\star$ & &$\star$& & N &[28]\\ 684{\tt bc95 }& & $\star$ & &$\star$& & N &[33]\\ 685{\tt cpw6b95 }& & $\star$ & &$\star$& & N &[36]\\ 686{\tt cpwb6k }& & $\star$ & &$\star$& & N &[36]\\ 687{\tt xm05 }& $\star$ & & &$\star$& & N &[37]\\ 688{\tt cm05 }& & $\star$ & &$\star$& & N &[37]\\ 689{\tt xm05-2x }& $\star$ & & &$\star$& & N &[38]\\ 690{\tt cm05-2x }& & $\star$ & &$\star$& & N &[38]\\ 691{\tt xctpssh }& & & &$\star$&$\star$& N &[29]\\ 692{\tt bb1k }& & & &$\star$&$\star$& N &[34]\\ 693{\tt mpw1b95 }& & & &$\star$&$\star$& N &[35]\\ 694{\tt mpwb1k }& & & &$\star$&$\star$& N &[35]\\ 695{\tt pw6b95 }& & & &$\star$&$\star$& N &[36]\\ 696{\tt pwb6k }& & & &$\star$&$\star$& N &[36]\\ 697{\tt m05 }& & & &$\star$&$\star$& N &[37]\\ 698{\tt vsxc }& & & &$\star$&$\star$& N &[39]\\ 699{\tt xvsxc }& $\star$ & & &$\star$& & N &[39]\\ 700{\tt cvsxc }& & $\star$ & &$\star$& & N &[39]\\ 701{\tt m06-L }& $\star$ & $\star$ & &$\star$&$\star$& N &[40]\\ 702{\tt xm06-L }& & & &$\star$& & N &[40]\\ 703{\tt cm06-L }& & $\star$ & &$\star$& & N &[40]\\ 704{\tt m06-hf }& & & &$\star$&$\star$& N &[41]\\ 705{\tt xm06-hf }& $\star$ & & &$\star$& & N &[41]\\ 706{\tt cm06-hf }& & $\star$ & &$\star$& & N &[41]\\ 707{\tt m06 }& & & &$\star$&$\star$& N &[42]\\ 708{\tt xm06 }& $\star$ & & &$\star$& & N &[42]\\ 709{\tt cm06 }& & $\star$ & &$\star$& & N &[42]\\ 710{\tt m06-2x }& & & &$\star$&$\star$& N &[42]\\ 711{\tt xm06-2x }& $\star$ & & &$\star$& & N &[42]\\ 712{\tt cm06-2x }& & $\star$ & &$\star$& & N &[42]\\ 713\hline 714{\tt beckehandh}& $\star$ & $\star$ & & & $\star$ & Y &[13]\\ 715{\tt b3lyp }& $\star$ & $\star$ & $\star$ & & $\star$ & Y &[14]\\ 716{\tt acm }& $\star$ & $\star$ & $\star$ & & $\star$ & Y &[14]\\ 717{\tt becke97 }& $\star$ & $\star$ & $\star$ & & $\star$ & N &[15]\\ 718{\tt becke97-1}& $\star$ & $\star$ & $\star$ & & $\star$ & N &[11]\\ 719{\tt becke97-2}& $\star$ & $\star$ & $\star$ & & $\star$ & N &[22]\\ 720{\tt becke97-3}& $\star$ & $\star$ & $\star$ & & $\star$ & N &[30]\\ 721{\tt becke97-d}& $\star$ & $\star$ & $\star$ & & $\star$ & N &[45]\\ 722{\tt becke98 }& $\star$ & $\star$ & $\star$ & & $\star$ & N &[16]\\ 723{\tt pbe0 }& $\star$ & $\star$ & $\star$ & & $\star$ & Y &[17]\\ 724{\tt mpw1k }& $\star$ & $\star$ & $\star$ & & $\star$ & Y &[25]\\ 725\hline 726\end{tabular} 727\caption{Table of available Exchange (X) and Correlation (C) functionals. 728GGA is the Generalized Gradient Approximation, and Meta refers to Meta-GGAs. 729The column {\em 2nd} refers to second derivatives of the 730energy with respect to nuclear position. } 731\label{tablexc} 732} 733\end{table} 734 735 736 737{\footnotesize 738\vspace{8.5cm} 739\begin{enumerate}\setlength{\itemsep}{-1\baselineskip} 740\setlength{\parsep}{-1\baselineskip} 741\item J.C.~Slater and K. H. Johnson, Phys. Rev. B {\bf 5}, 844 (1972)\\ 742\item S.J.~Vosko, L.~Wilk and M.~Nusair, 743Can.~J.~Phys.~{\bf 58}, 1200 (1980)\\ 744\item J.~P.~Perdew and A.~Zunger, Phys.~Rev.~B {\bf23}, 5048 745(1981). \\ 746\item J.P.~Perdew and Y.~Wang, Phys. Rev. B {\bf 45}, 74713244 (1992)\\ 748\item A.D.~Becke, Phys. Rev. A {\bf 88}, 3098 (1988)\\ 749\item J.P. Perdew, J.A.~Chevary, S.H.~Vosko, K.A.~Jackson, 750M.R.~Pederson, D.J.~Singh and C.~Fiolhais, Phys. Rev. B {\bf 46}, 6671 751(1992)\\ 752\item J.P.~Perdew, K.~Burke and M.~Ernzerhof, 753Phys. Rev. Lett. {\bf 77}, 3865 (1996); {\bf 78 }, 1396 (1997)\\ 754\item P.W.Gill , Mol. Phys. {\bf 89}, 433 (1996)\\ 755\item J.~P.~Perdew, Phys.~Rev.~B {\bf33}, 8822 (1986)\\ 756\item C.~Lee, W.~Yang and R.~G.~Parr, Phys.~Rev.~B {\bf 37}, 785 757(1988)\\ 758\item F.A.Hamprecht, A.J.Cohen, D.J.Tozer and N.C\\ Handy, 759 J. Chem. Phys. {\bf 109}, 6264 (1998)\\ 760\item A.D.Boese, N.L.Doltsinis, N.C.Handy and 761M.Sprik. J. Chem. Phys. {\bf 112}, 1670 (2000)\\ 762\item A.D.~Becke, J.~Chem.~Phys. {\bf 98}, 1372 (1992)\\ 763\item A.D.~Becke, J.~Chem.~Phys.~{\bf 98}, 5648 (1993)\\ 764\item A.D.Becke, J. Chem. Phys. {\bf 107}, 8554 (1997)\\ 765\item H.L.Schmider and A.D.~Becke, J.~Chem.~Phys.~{\bf 108}, 7669624 (1998)\\ 767\item C.Adamo and V.Barone, J.~Chem.~Phys. {\bf 110}, 6158 (1998)\\ 768\item A.J.Cohen and N.C. Handy, Chem. Phys. Lett. {\bf 316}, 160 (2000)\\ 769\item A.D.Boese, N.C.Handy, J. Chem. Phys. {\bf 114}, 5497 770(2001)\\ 771\item N.C.Handy, A.J. Cohen, Mol. Phys. {\bf 99}, 403 (2001)\\ 772\item G. Menconi, P.J. Wilson, D.J. Tozer, 773J. Chem. Phys {\bf 114}, 3958 (2001)\\ 774\item P.J. Wilson, T.J. Bradley, D.J. Tozer, J. Chem. Phys {\bf 115}, 7759233 (2001)\\ 776\item C. Adamo and V. Barone, J.~Chem.~Phys. {\bf 108}, 664 (1998); Y. Zhao and D.G. Truhlar, J. Phys. Chem. A {\bf 109}, 5656 (2005)\\ 777\item M.Filatov and W.Thiel, Mol.Phys. {\bf 91}, 847 (1997). 778 M.Filatov and W.Thiel, Int.J.Quantum Chem. {\bf 62}, 603 (1997)\\ 779\item B.J. Lynch, P.L. Fast, M. Harris and D.G. Truhlar, J. Phys. Chem. A 780{\bf 104}, 4811(2000)\\ 781\item J.P.~Perdew, S.~Kurth, A.~Zupan and P.~Blaha, 782Phys. Rev. Lett. {\bf 82}, 2544 (1999)\\ 783\item A.~D.~Boese, A.~Chandra, J.~M.~L.~Martin and D.~Marx, 784 J. Chem. Phys. {\bf 119}, 5965 (2003)\\ 785\item J. Tao,J.Perdew, V. Staroverov and G. Scuseria, 786Phys. Rev. Let. {\bf 91}, 146401-1 (2003)\\ 787\item V. Staroverov, G. Scuseria, 788J. Tao and J.Perdew, J. Chem.Phys. {\bf 119}, 12129 (2003)\\ 789\item T.W. Keal, D.J. Tozer, J. Chem. Phys {\bf 123}, 790121103 (2005)\\ 791\item T. Tsuneda, T. Suzumura and K. Hirao, 792J. Chem Phys. {\bf 110}, 10664 (1999)\\ 793\item T. Tsuneda, T. Suzumura and K. Hirao, 794J. Chem Phys. {\bf 111}, 5656 (1999) \\ 795\item B. Hammer, L. B. Hansen and J. N{\o}rskov , Phys. Rev. B {\bf 58}, 7413 (1999)\\ 796\item Y. Zhang and W. Yang, Phys. Rev. Letters {\bf 80}, 890 (1998)\\ 797\item A. D. Becke, J. Chem. Phys. {\bf 104}, 1040 (1996)\\ 798\item Y. Zhao and D. G. Truhlar, J. Phys. Chem. A {\bf 108}, 2715 (2004)\\ 799\item Y. Zhao and D. G. Truhlar, J. Phys. Chem. A {\bf 108}, 6908 (2004)\\ 800\item Y. Zhao and D. G. Truhlar, J. Phys. Chem. A {\bf 109}, 5656 (2005)\\ 801\item Y. Zhao, N. E. Schultz and D. G. Truhlar, J. Chem. Phys. {\bf 123}, 161103 (2005)\\ 802\item Y. Zhao, N. E. Schultz and D. G. Truhlar, J. Chem. Theory Comput. {\bf 2}, 364 (2006)\\ 803\item T. Van Voorhis, G. E. Scuseria, J. Chem. Phys. {\bf 109}, 400 (1998)\\ 804\item Y. Zhao, D. G. Truhlar, J. Chem. Phys. {\bf 125}, 194101 (2006)\\ 805\item Y. Zhao, D. G. Truhlar, J. Phys. Chem. A. {\bf 110}, 13126 (2006)\\ 806\item Y. Zhao, D. G. Truhlar, Theor. Chem. Acc. (2006)\\ 807\item S. Grimme, J. Comp. Chem. \bf{27}, 1787 (2006). 808 809\end{enumerate} 810} 811%\end{table} 812\onecolumn 813 814\section{{\tt ITERATIONS} --- Number of SCF iterations} 815 816\begin{verbatim} 817 ITERATIONS <integer iterations default 30> 818\end{verbatim} 819 820The default optimization in the DFT module is to iterate on the 821Kohn-Sham (SCF) equations for a specified number of iterations 822(default 30). The keyword that controls this optimization 823is \verb+ITERATIONS+, and has the following general form, 824 825\begin{verbatim} 826 iterations <integer iterations default 30> 827\end{verbatim} 828 829The optimization procedure will stop when the specified number of 830iterations is reached or convergence is met. See an example 831that uses this directive in section \ref{sec:DFTsample}. 832 833\section{{\tt CONVERGENCE} --- SCF Convergence Control} 834\label{sec:dftconv} 835 836\begin{verbatim} 837 CONVERGENCE [energy <real energy default 1e-6>] \ 838 [density <real density default 1e-5>] \ 839 [gradient <real gradient default 5e-4>] \ 840 [hl_tol <real hl_tol default 0.1>] 841 [dampon <real dampon default 0.0>] \ 842 [dampoff <real dampoff default 0.0>] \ 843 [ncydp <integer ncydp default 2>] \ 844 [ncyds <integer ncyds default 30>] \ 845 [ncysh <integer ncysh default 30>] \ 846 [damp <integer ndamp default 0>] [nodamping] \ 847 [diison <real diison default 0.0>] \ 848 [diisoff <real diisoff default 0.0>] \ 849 [(diis [nfock <integer nfock default 10>]) || nodiis] \ 850 [levlon <real levlon default 0.0>] \ 851 [levloff <real levloff default 0.0>] \ 852 [(lshift <real lshift default 0.5>) || nolevelshifting] \ 853 [rabuck [n_rabuck <integer n_rabuck default 25>]] 854\end{verbatim} 855 856Convergence is satisfied by meeting any or all of three criteria; 857\begin{itemize} 858\item convergence of the total energy; this is defined to be when the 859 total DFT energy at iteration N and at iteration N-1 differ by a value less 860 than some value (the default is 1e-6). This value can be modified 861 using the key word, 862\begin{verbatim} 863 CONVERGENCE energy <real energy default 1e-6> 864\end{verbatim} 865 866\item convergence of the total density; this is defined to be when the 867 total DFT density matrix at iteration N and at iteration N-1 have a 868 RMS difference less than some value (the default is 1e-5). This value can be modified 869 using the key word, 870\begin{verbatim} 871 CONVERGENCE density <real density default 1e-5> 872\end{verbatim} 873 874\item convergence of the orbital gradient; this is defined to be when the 875 DIIS error vector becomes less than some value (the default is 876 5e-4). This value can be modified using the key word, 877\begin{verbatim} 878 CONVERGENCE gradient <real gradient default 5e-4> 879\end{verbatim} 880\end{itemize} 881 882The default optimization strategy is to immediately begin direct 883inversion of the iterative subspace\footnote {P.~Pulay, Chem.\ Phys.\ 884 Lett.\ {\bf 73}, 393 (1980) and P.~Pulay, J.~Comp.~Chem.~{\bf 3}, 885 566 (1982)}. Damping is also initiated (using 70\% of the previous 886density) for the first 2 iteration. In addition, if the HOMO - LUMO 887gap is small and the Fock matrix somewhat diagonally dominant, then 888level-shifting is automatically initiated. There are a variety of ways 889to customize this procedure to whatever is desired. 890 891An alternative optimization strategy is to specify, by using the change 892in total energy (from iterations when N and N-1), when to turn 893damping, level-shifting, and/or DIIS on/off. Start and stop keywords for 894each of these is available as, 895\begin{verbatim} 896 CONVERGENCE [dampon <real dampon default 0.0>] \ 897 [dampoff <real dampoff default 0.0>] \ 898 [diison <real diison default 0.0>] \ 899 [diisoff <real diisoff default 0.0>] \ 900 [levlon <real levlon default 0.0>] \ 901 [levloff <real levloff default 0.0>] 902\end{verbatim} 903 904So, for example, damping, DIIS, and/or level-shifting can be turned 905on/off as desired. 906 907Another strategy can be to simply specify how many iterations (cycles) you wish 908each type of procedure to be used. The necessary keywords to control 909the number of damping cycles (ncydp), the number of DIIS cycles 910(ncyds), and the number of level-shifting cycles (ncysh) are input as, 911\begin{verbatim} 912 CONVERGENCE [ncydp <integer ncydp default 2>] \ 913 [ncyds <integer ncyds default 30>] \ 914 [ncysh <integer ncysh default 0>] 915\end{verbatim} 916 917The amount of damping, level-shifting, time at which level-shifting is 918automatically imposed, and Fock matrices used in the DIIS 919extrapolation can be modified by the following keywords 920\begin{verbatim} 921 CONVERGENCE [damp <integer ndamp default 0>] \ 922 [diis [nfock <integer nfock default 10>]] \ 923 [lshift <real lshift default 0.5>] \ 924 [hl_tol <real hl_tol default 0.1>]] 925\end{verbatim} 926 927Damping is defined to be the percentage of the previous iterations 928density mixed with the current iterations density. So, for example 929\begin{verbatim} 930 CONVERGENCE damp 70 931\end{verbatim} 932would mix 30\% of the current iteration density with 70\% of the 933previous iteration density. 934 935Level-Shifting\footnote {M.F.~Guest and 936V.R.~Saunders, Mol.~Phys.~{\bf 28}, 819 (1974)} is defined as the 937amount of shift applied to the diagonal elements of the unoccupied 938block of the Fock matrix. The shift is specified by the 939keyword \verb+lshift+. For example the directive, 940\begin{verbatim} 941 CONVERGENCE lshift 0.5 942\end{verbatim} 943causes the diagonal elements of the Fock matrix 944corresponding to the virtual orbitals to be shifted by 0.5 a.u. 945By default, this level-shifting procedure is switched on whenever the 946HOMO-LUMO gap is small. Small is defined by default to be 0.05 au but 947can be modified by the directive \verb+hl_tol+. An example of 948changing the HOMO-LUMO gap tolerance to 0.01 would be, 949\begin{verbatim} 950 CONVERGENCE hl_tol 0.01 951\end{verbatim} 952 953Direct inversion of the iterative subspace with extrapolation of up to 95410 Fock matrices is a default optimization procedure. For large 955molecular systems the amount of available memory may preclude the ability to 956store this number of N**2 arrays in global memory. The user may then 957specify the number of Fock matrices to be used in the extrapolation 958(must be greater than three (3) to be effective). To set the number of 959Fock matrices stored and used in the extrapolation procedure to 3 960would take the form, 961\begin{verbatim} 962 CONVERGENCE diis 3 963\end{verbatim} 964 965The user has the ability to simply turn off any optimization 966procedures deemed undesirable with the obvious keywords, 967\begin{verbatim} 968 CONVERGENCE [nodamping] [nodiis] [nolevelshifting] 969\end{verbatim} 970 971 972For systems where the initial guess is very poor, the user can try the 973method described in 974\footnote{A.~D.~Rabuck and G.~E.~Scuseria, J.~Chem.~Phys {\bf 110},695 975(1999)} 976that makes use of {\bf fractional occupation} of the orbital levels during 977the initial cycles of the SCF convergence. The input has the following form 978 979\begin{verbatim} 980 CONVERGENCE rabuck [n_rabuck <integer n_rabuck default 25>]] 981\end{verbatim} 982 983where the optional value {\tt n\_rabuck} determines the number of SCF 984cycles during which the method will be active. For example, to 985set equal to 30 the number of cycles where the Rabuck method is 986active, you need to use the following line 987\begin{verbatim} 988 CONVERGENCE rabuck 30 989\end{verbatim} 990 991\section{{\tt CDFT} --- Constrained DFT} 992\label{cdft} 993 994This option enables the constrained DFT formalism by Wu and Van Voorhis described 995in the paper: Q. Wu, T. Van Voorhis, Phys. Rev. A {\bf 72}, 024502 (2005). 996 997\begin{verbatim} 998 CDFT <integer fatom1 latom1> [<integer fatom2 latom2>] (charge||spin <real constaint_value>) \ 999 [pop (becke||mulliken||lowdin) default lowdin] 1000\end{verbatim} 1001 1002Variables fatom1 and latom1 define the first and last atom of the group of atoms to which 1003the constaint will be applied. Therefore the atoms in the same group should be placed 1004continuously in the geometry input. If fatom2 and latom2 are specified, the difference between 1005group 1 and 2 (i.e. 1-2) is constrained. 1006 1007The constraint can be either on the charge or the spin density (\# of alpha - beta electrons) with 1008a user specified constaint\_value. Note: No gradients have been implemented for the spin constaints 1009case. Geometry optimizations can only be performed using the charge constaint. 1010 1011To calculate the charge or spin density, the Becke, Mulliken, and Lowdin population schemes can be 1012used. The Lowdin scheme is default while the Mulliken scheme is not recommended. If basis sets with 1013many diffuse functions are used, the Becke population scheme is recommended. 1014 1015Multiple constaints can be defined simultaniously by defining multiple {\tt cdft} lines in the input. 1016The same population scheme will be used for all constaints and only needs to be specified once. If 1017multiple population options are defined, the last one will be used. When there are convergence 1018problems with multiple constaints, the user is advised to do one constraint first and to use the 1019resulting orbitals for the next step of the constained calculations. 1020 1021It is best to put "convergence nolevelshifting" in the dft directive to avoid issues with gradient 1022calculations and convergence in CDFT. Use orbital swap to get a broken-symmetry solution. 1023 1024An input example is given below. 1025 1026\begin{verbatim} 1027geometry 1028symmetry 1029C 0.0 0.0 0.0 1030O 1.2 0.0 0.0 1031C 0.0 0.0 2.0 1032O 1.2 0.0 2.0 1033end 1034 1035basis 1036* library 6-31G* 1037end 1038 1039dft 1040 xc b3lyp 1041 convergence nolevelshifting 1042 odft 1043 mult 1 1044 vectors swap beta 14 15 1045 cdft 1 2 charge 1.0 1046end 1047task dft 1048\end{verbatim} 1049 1050\section{{\tt SMEAR} --- Fractional Occupation of the Molecular Orbitals} 1051\label{smear} 1052 1053The {\tt \bf SMEAR} keyword is useful in cases with many degenerate states 1054near the HOMO (eg metallic clusters) 1055 1056\begin{verbatim} 1057 SMEAR <real smear default 0.001> 1058\end{verbatim} 1059 1060This option allows fractional occupation of the molecular orbitals. 1061A Gaussian broadening function of exponent {\tt smear} is used as described in 1062the paper: 1063R.W. Warren and B.I. Dunlap, Chem. Phys. Letters {\bf 262}, 384 (1996).\\ 1064The user must be aware that an additional energy term is added to the total 1065energy in order to have 1066energies and gradients consistent. 1067 1068 1069\section{{\tt GRID} --- Numerical Integration of the XC Potential} 1070\label{grgrid} 1071\begin{verbatim} 1072 GRID [(xcoarse||coarse||medium||fine||xfine) default medium] \ 1073 [(gausleg||lebedev ) default lebedev ] \ 1074 [(becke||erf1||erf2||ssf) default erf1] \ 1075 [(euler||mura||treutler) default mura] \ 1076 [rm <real rm default 2.0>] \ 1077 [nodisk] 1078\end{verbatim} 1079 1080A numerical integration is necessary for the evaluation of the 1081exchange-correlation contribution to the density functional. The 1082default quadrature used for the numerical integration is an 1083Euler-MacLaurin scheme for the radial components (with a modified 1084Mura-Knowles transformation) 1085and a Lebedev 1086scheme for the angular components. Within this numerical 1087integration procedure various levels of accuracy have been defined and 1088are available to the user. The user can specify the level of accuracy 1089with the keywords; xcoarse, coarse, medium, fine, and xfine. The 1090default is medium. 1091 1092\begin{verbatim} 1093 GRID [xcoarse||coarse||medium||fine||xfine] 1094\end{verbatim} 1095 1096Our intent is to have a numerical integration scheme which would give 1097us approximately the accuracy defined below regardless of molecular 1098composition. 1099\begin{center} 1100 \begin{tabular}[right]{|c|c|} \hline 1101Keyword & {\tt Total Energy Target Accuracy} \\ \hline 1102{\tt xcoarse} & $1x10^{-4}$ \\ \hline 1103{\tt coarse} & $1x10^{-5}$ \\ \hline 1104{\tt medium} & $1x10^{-6}$ \\ \hline 1105{\tt fine} & $1x10^{-7}$ \\ \hline 1106{\tt xfine} & $1x10^{-8}$ \\ \hline 1107 \end{tabular} \\ 1108\end{center} 1109 1110In order to determine the level of radial and angular quadrature needed 1111to give us the target accuracy we computed total DFT energies 1112at the LDA level of theory for many 1113homonuclear atomic, diatomic and triatomic systems in rows 1-4 of the 1114periodic table. In each case all bond lengths were set to twice the 1115Bragg-Slater radius. The total DFT energy of the system was computed 1116using the converged SCF density with atoms having radial shells 1117ranging from 35-235 (at fixed 48/96 angular quadratures) and angular 1118quadratures of 12/24-48/96 (at fixed 235 radial shells). The error of 1119the numerical integration was determined by comparison to a ``best'' 1120or most accurate calculation in which a grid of 235 radial points 48 1121theta and 96 phi angular points on each atom was used. This 1122corresponds to approximately 1 million points per atom. The following 1123tables were empirically determined to give the desired target accuracy 1124for DFT total energies. These tables below show the number of radial and 1125angular shells which the DFT module will use for for a given atom 1126depending on the row it is in (in the periodic table) and the desired 1127accuracy. Note, differing atom types in a given molecular system will 1128most likely have differing associated numerical grids. The intent is 1129to generate the desired energy accuracy (with utter disregard for speed). 1130 1131%{\bf Important note to users.} We clearly understand that the default 1132%(Euler-MacLaurin/Gauss-Legendre) grids are large and result in slow 1133%%construction of the numerical components of the Kohn-Sham equations. 1134%%Alternatively, we have provided access to two-dimensional Lebedev 1135%%angular quadratures which can be used in many cases to substantially 1136%%reduce the number of grid points per atom while keeping good accuracy. 1137%%We have not yet had the opportunity to benchmark the Lebedev angular 1138%%quadratures to the same extent that we have for the Gauss-Legendre. 1139%%We therefore do not have default Lebedev quadratures for specific 1140%%target accuracy for all elements of the periodic table. If the user 1141%%wants to significantly decrease CPU time to solution it is suggested 1142%%that a few prototype benchmark calculations be done using various 1143%%Lebedev quadratures (which we describe below) while monitoring the 1144%%numerically integrated density and total energies for the molecular 1145%%systems of interest. For many examples we have observed speed-ups of 1146%%two or more for the same numerical accuracy when using Lebedev rather 1147%%than the default Gauss-Legendre quadrature. In addition, we have 1148%%observed that with Lebedev angular quadratures a reduction in the 1149%%number of radial shells (perhaps by as much as 30\%) might be possible 1150%%while continuing to provide the same level of accuracy. 1151%%% g94/Sg-1 grid ssf lebedev gausleg 50 8 1152%%% g98/fine grid ssf lebedev gausleg 75 11 1153 1154 1155\begin{table}[h] 1156\begin{center} 1157\caption{Program default number of radial and angular shells empirically determined for Row 1 atoms 1158 (Li $\rightarrow$ F) to reach the desired accuracies.} 1159 1160\vspace{.2in} 1161 1162 \begin{tabular}[right]{|c|c|c|} \hline 1163Keyword & {\tt Radial} & {\tt Angular} \\ \hline 1164{\tt xcoarse} & 21 & 194 \\ \hline 1165{\tt coarse} & 35 & 302 \\ \hline 1166{\tt medium} & 49 & 434 \\ \hline 1167{\tt fine} & 70 & 590 \\ \hline 1168{\tt xfine} &100 &1202 \\ \hline 1169 \end{tabular} \\ 1170\end{center} 1171\end{table} 1172 1173\begin{table}[h] 1174\begin{center} 1175\caption{Program default number of radial and angular shells empirically determined for Row 2 atoms 1176 (Na $\rightarrow$ Cl) to reach the desired accuracies.} 1177 1178\vspace{.2in} 1179 1180 \begin{tabular}[right]{|c|c|c|} \hline 1181Keyword & {\tt Radial} & {\tt Angular} \\ \hline 1182{\tt xcoarse} & 42 & 194 \\ \hline 1183{\tt coarse} & 70 & 302 \\ \hline 1184{\tt medium} & 88 & 434 \\ \hline 1185{\tt fine} &123 & 770 \\ \hline 1186{\tt xfine} &125 &1454 \\ \hline 1187 \end{tabular} \\ 1188\end{center} 1189\end{table} 1190 1191\begin{table}[h] 1192\begin{center} 1193\caption{Program default number of radial and angular shells empirically determined for Row 3 atoms 1194 (K $\rightarrow$ Br) to reach the desired accuracies.} 1195 1196\vspace{.2in} 1197 1198 \begin{tabular}[right]{|c|c|c|} \hline 1199Keyword & {\tt Radial} & {\tt Angular} \\ \hline 1200{\tt xcoarse} & 75 & 194 \\ \hline 1201{\tt coarse} & 95 & 302 \\ \hline 1202{\tt medium} &112 & 590 \\ \hline 1203{\tt fine} &130 & 974 \\ \hline 1204{\tt xfine} &160 &1454 \\ \hline 1205 \end{tabular} \\ 1206\end{center} 1207\end{table} 1208 1209\begin{table}[h] 1210\begin{center} 1211\caption{Program default number of radial and angular shells empirically determined for Row 4 atoms 1212 (Rb $\rightarrow$ I) to reach the desired accuracies.} 1213 1214\vspace{.2in} 1215 1216 \begin{tabular}[right]{|c|c|c|} \hline 1217Keyword & {\tt Radial} & {\tt Angular} \\ \hline 1218{\tt xcoarse} & 84 &194 \\ \hline 1219{\tt coarse} &104 &302 \\ \hline 1220{\tt medium} &123 &590 \\ \hline 1221{\tt fine} &141 &974 \\ \hline 1222{\tt xfine} &205 &1454 \\ \hline 1223 \end{tabular} \\ 1224\end{center} 1225\end{table} 1226 1227\clearpage 1228 1229\subsection{Angular grids} 1230 1231In addition to the simple keyword specifying the desired accuracy as 1232described above, the user has the option of specifying a custom 1233quadrature of this type in which ALL atoms have the same grid 1234specification. This is accomplished by using the \verb+gausleg+ keyword. 1235 1236\paragraph{Gauss-Legendre angular grid} 1237 1238\begin{verbatim} 1239 GRID gausleg <integer nradpts default 50> <integer nagrid default 10> 1240\end{verbatim} 1241 1242In this type of grid, the number of phi points is twice the number of 1243theta points. So, for example, a specification of, 1244\begin{verbatim} 1245 GRID gausleg 80 20 1246\end{verbatim} 1247would be interpreted as 80 radial points, 20 theta points, and 40 1248phi points per center (or 64000 points per center before pruning). 1249 1250\paragraph{Lebedev angular grid} 1251 1252A second quadrature is the Lebedev 1253scheme for the angular components\footnote{The subroutine 1254for the Lebedev grid was derived from a routine supplied by M.~Caus\`a 1255of the University of Torino and from the grid points supplied by 1256D.N.~Laikov from Moscow State University.}. 1257Within this numerical integration procedure various levels 1258of accuracy have also been defined and are available to the user. 1259The input for this type of grid takes the form, 1260\begin{verbatim} 1261 GRID lebedev <integer radpts > <integer iangquad > 1262\end{verbatim} 1263In this context the variable iangquad specifies a certain number of 1264angular points as indicated by the table below.\footnote{ 1265V.I. Lebedev and D.N. Laikov, Doklady Mathematics {\bf 366}, 741 1266(1999). 1267} 1268\begin{table}[htp] 1269\begin{center} 1270\begin{tabular}[right]{|c|rr|} \hline 1271$IANGQUAD$ & $N_{angular}$ & $l$\\ 1272\hline 1273 1& 38& 9 \\ 1274 2& 50& 11 \\ 1275 3& 74 & 13 \\ 1276 4& 86 & 15 \\ 1277 5& 110 & 17 \\ 1278 6& 146 & 19 \\ 1279 7& 170 & 21 \\ 1280 8& 194 & 23 \\ 1281 9& 230 & 25 \\ 128210& 266& 27 \\ 128311& 302& 29 \\ 128412& 350& 31 \\ 128513& 434& 35 \\ 128614& 590& 41 \\ 128715& 770& 47 \\ 128816& 974& 53 \\ 128917& 1202& 59 \\ 129018& 1454& 65 \\ 129119& 1730& 71 \\ 129220& 2030& 77 \\ 129321& 2354& 83 \\ 129422& 2702& 89 \\ 129523& 3074& 95 \\ 129624& 3470& 101 \\ 129725& 3890& 107 \\ 129826& 4334& 113 \\ 129927& 4802& 119 \\ 130028& 5294& 125 \\ 130129& 5810& 131 \\ 1302\hline 1303\end{tabular} 1304\end{center} 1305\caption{List of Lebedev quadratures} 1306\end{table} 1307Therefore the user can specify any number of radial points along with 1308the level of angular quadrature (1-29). 1309 1310The user can also specify grid parameters specific for a given atom type: 1311parameters that must be supplied are: atom tag and number of radial points. 1312As an example, here is a grid input line for the water molecule 1313\begin{verbatim} 1314grid lebedev 80 11 H 70 8 O 90 11 1315\end{verbatim} 1316 1317 1318% Delley weights do not work .. need to replace with SSWs after input is 1319%changed 1320%The user also has the option of choosing one of two types of spatial weights 1321%implemented in the numerical integration of the XC terms; Delley and Becke. 1322 1323\clearpage 1324\subsection{Partitioning functions} 1325 1326\begin{verbatim} 1327 GRID [(becke||erf1||erf2||ssf) default erf1] 1328\end{verbatim} 1329 1330 1331\begin{description} 1332\item[\tt becke] A. D. Becke, J. Chem. Phys. {\bf 88}, 1053 (1988). 1333\item[\tt ssf] R.E.Stratmann, G.Scuseria and M.J.Frisch, 1334Chem. Phys. Lett. {\bf 257}, 213 (1996). 1335\item[\tt erf1] modified ssf 1336\item[\tt erf2] modified ssf 1337\end{description} 1338 1339Erf$n$ partioning functions 1340 1341\begin{eqnarray*} 1342 w_A(r) & = & \prod_{B\neq A}\frac{1}{2} \left[1 \ - \ 1343erf(\mu^\prime_{AB})\right] \\ 1344 \mu^\prime_{AB} & = & \frac{1}{\alpha} \ \frac{\mu_{AB}}{(1-\mu_{AB}^2)^n}\\ 1345 \mu_{AB} & = & \frac{{\mathbf r}_A - {\mathbf r}_B} 1346{\left|{\mathbf r}_A - {\mathbf r}_B \right|} 1347\end{eqnarray*} 1348 1349 1350 1351 1352\subsection{Radial grids} 1353 1354\begin{verbatim} 1355 GRID [[euler||mura||treutler] default mura] 1356\end{verbatim} 1357 1358\begin{description} 1359\item[\tt euler] Euler-McLaurin quadrature wih the transformation 1360 devised by 1361C.W. Murray, N.C. Handy, and G.L. Laming, 1362Mol. Phys.{\bf 78}, 997 (1993). 1363 \\ 1364\item[\tt mura] Modification of the Murray-Handy-Laming scheme by 1365M.E.Mura and P.J.Knowles, J Chem Phys {\bf 104}, 9848 1366(1996) (we are not using the scaling factors proposed 1367in this paper).\\ 1368\item[\tt treutler] Gauss-Chebyshev using the transformation suggested 1369 by O.Treutler and R.Alrhichs, J.Chem.Phys {\bf 102}, 346 (1995).\\ 1370\end{description} 1371 1372\subsection{Disk usage for Grid} 1373 1374\begin{verbatim} 1375 NODISK 1376\end{verbatim} 1377 1378This keyword turns off storage of grid points and weights on disk. 1379 1380 1381\section{{\tt TOLERANCES} --- Screening tolerances} 1382 1383\begin{verbatim} 1384 TOLERANCES [[tight] [tol_rho <real tol_rho default 1e-10>] \ 1385 [accCoul <integer accCoul default 8>] \ 1386 [radius <real radius default 25.0>]] 1387\end{verbatim} 1388% [accQrad <integer accQrad default 12>] \ 1389% [accAOfunc <integer accAOfunc default 20>] \ 1390% [accXCfunc <integer accXCfunc default 20>] \ 1391% [accCDfunc <integer accAOfunc default 20>] \ 1392% 1393%{\bf JEFF: tight needs explanation} 1394% 1395The user has the option of controlling screening for the tolerances in 1396the integral evaluations for the DFT module. In most applications, 1397the default values will be adequate for the calculation, but different 1398values can be specified in the input for the DFT module using the 1399keywords described below. 1400 1401%The input to define a screening tolerance for evaluation of the AO 1402%Gaussian functions is specified with the keyword \verb+accAOfunc+, as 1403%follows, 1404%\begin{verbatim} 1405% TOLERANCES accAOfunc <integer accAOfunc> 1406%\end{verbatim} 1407%A Gaussian Type Function of the Orbital basis set is 1408%evaluated at a point $r_i$ if its value it is greater than 1409%$exp(-${\tt accAOfunc}$)$ ; 1410%the default value is set to $-ln(\Delta E)$, where $\Delta E$ is the desired 1411%accuracy on energy. 1412%A Gaussian orbital basis (AO) function with exponent $\zeta$ 1413%and radial factor $e^{-\zeta\cdot r_i^2}$ is 1414%evaluated at a point $r_i$ only if 1415%$\zeta\cdot r_i^2$ is less than the value specified for ${\tt accAOfunc}$. 1416 1417%The input to define a screening tolerance for evaluation of the exchange- 1418%correlation (XC) Gaussian fitting functions is specified with the 1419%keyword \verb+accXCfunc+, as follows, 1420%\begin{verbatim} 1421% TOLERANCES accXCfunc <integer accXCfunc default 20> 1422%\end{verbatim} 1423%An exchange-correlation (XC) fitting function with exponent $\zeta$ 1424%and radial factor $e^{-\zeta\cdot r_i^2}$ is 1425%evaluated at a point $r_i$ only if 1426%$\zeta\cdot r_i^2$ is less than the value specified for ${\tt accXCfunc}$. 1427% 1428%The input to define a screening tolerance for evaluation of the 1429%charge-density (CD) Gaussian fitting functions is specified with the 1430%keyword \verb+accCDfunc+, as follows, 1431%\begin{verbatim} 1432% TOLERANCES accCDfunc <integer accCDfunc default 20> 1433%\end{verbatim} 1434%A charge-density (CD) fitting function with exponent $\zeta$ 1435%and radial factor $e^{-\zeta\cdot r_i^2}$ is evaluated at a 1436%point $r_i$ only if $\zeta\cdot r_i^2$ is less than the value 1437%specified for ${\tt accCDfunc}$. 1438 1439The input 1440parameter {\tt accCoul} is used to define the tolerance in Schwarz 1441screening for the Coulomb integrals. Only integrals with estimated 1442values greater than $10^{(-{\tt accCoul})}$ are evaluated. 1443 1444\begin{verbatim} 1445 TOLERANCES accCoul <integer accCoul default 8> 1446\end{verbatim} 1447 1448%The user also has the option of specifying the radial quadrature 1449%grid cut-off for the DFT calculation, using the keyword 1450%\verb+accQrad+. The input line for this option is as follows, 1451%\begin{verbatim} 1452% TOLERANCES accQrad <integer accQrad default 12> 1453%\end{verbatim} 1454 1455%The value entered for \verb+accQrad+ is the cutoff distance, in bohr, for grid 1456%points around a given center or atom. Grid points that lie more than 1457%\verb+accQrad+ bohr from the center or atom are neglected. 1458 1459Screening away needless computation of the XC functional (on the grid) 1460due to negligible density is also possible with the use of, 1461\begin{verbatim} 1462 TOLERANCES tol_rho <real tol_rho default 1e-10> 1463\end{verbatim} 1464XC functional computation is bypassed if the corresponding density 1465elements are less than \verb+tol_rho+. 1466 1467A screening parameter, \verb+radius+, used in the screening of the 1468Becke or Delley spatial weights is also available as, 1469\begin{verbatim} 1470 TOLERANCES radius <real radius default 25.0> 1471\end{verbatim} 1472where radius is the cutoff value in bohr. 1473 1474The tolerances as discussed previously are insured at convergence. 1475More sleazy tolerances are invoked early in the iterative process 1476which can speed things up a bit. This can also be problematic at 1477times because it introduces a discontinuity in the convergence 1478process. To avoid use of initial sleazy tolerances the user can 1479invoke the \verb+tight+ option: 1480 1481\begin{verbatim} 1482 TOLERANCES tight 1483\end{verbatim} 1484 1485This option sets all tolerances to their 1486default/user specified values at the very first iteration. 1487 1488 1489\section{{\tt DIRECT}, {\tt SEMIDIRECT} and {\tt NOIO} --- Hardware Resource Control} 1490\begin{verbatim} 1491 DIRECT||INCORE 1492 SEMIDIRECT [filesize <integer filesize default disksize>] 1493 [memsize <integer memsize default available>] 1494 [filename <string filename default $file_prefix.aoints$>] 1495 NOIO 1496\end{verbatim} 1497 1498\sloppy 1499 1500The inverted charge-density and exchange-correlation matrices 1501for a DFT calculation are normally written to disk storage. The user 1502can prevent this by specifying the keyword \verb+noio+ within the 1503input for the DFT directive. The input to exercise this option is 1504as follows, 1505\begin{verbatim} 1506 noio 1507\end{verbatim} 1508If this keyword is encountered, then the two matrices (inverted 1509charge-density and exchange-correlation) are computed ``on-the-fly'' 1510whenever needed. 1511 1512The \verb+INCORE+ option is always assumed to be true but can be 1513overridden with the option \verb+DIRECT+ in which case all integrals 1514are computed ``on-the-fly''. 1515 1516The \verb+SEMIDIRECT+ option controls caching of integrals. A full 1517description of this option is described in User Manual~\ref{sec:max}. 1518Some functionality which is only compatible with the \verb+DIRECT+ 1519option will not, at present, work when using \verb+SEMIDIRECT+. 1520 1521\fussy 1522 1523\section{{\tt ODFT} and {\tt MULT} --- Open shell systems} 1524\begin{verbatim} 1525 ODFT 1526 MULT <integer mult default 1> 1527\end{verbatim} 1528 1529Both {\sl closed-shell} and {\sl open-shell} systems can be studied using 1530the DFT module. Specifying the keyword \verb+MULT+ within the \verb+DFT+ 1531directive allows the user to define the spin multiplicity of the system. 1532The form of the input line is as follows; 1533\begin{verbatim} 1534 MULT <integer mult default 1> 1535\end{verbatim} 1536When the keyword \verb+MULT+ is specified, the user can define the integer 1537variable \verb+mult+, where \verb+mult+ is equal to the number of alpha 1538electrons minus beta electrons, plus 1. 1539 1540The keyword \verb+ODFT+ is unnecessary except in the context 1541of forcing a singlet system to be computed as an open shell 1542system (i.e., using a spin-unrestricted wavefunction). 1543 1544\section{{\tt SIC} --- Self-Interaction Correction} 1545 1546\begin{verbatim} 1547sic [perturbative || oep || oep-loc <default perturbative>] 1548\end{verbatim} 1549 1550The Perdew and Zunger (see J. P. Perdew and A. Zunger, Phys. Rev. B 23, 15515048 (1981)) method to remove the self-interaction contained in many 1552exchange-correlation functionals has been implemented with the 1553Optimized Effective Potential method 1554(see R. T. Sharp and G. K. Horton, Phys. Rev. {\bf 90}, 317 (1953), 1555J. D. Talman and W. F. Shadwick, Phys. Rev. A {\bf 14}, 36 (1976)) 1556within the Krieger-Li-Iafrate approximation (J. B. Krieger, Y. Li, 1557and G. J. Iafrate, Phys. Rev. A {\bf 45}, 101 (1992); {\bf 46}, 5453 (1992); 155847, 165 (1993)) 1559Three variants of these methods are included in NWChem: 1560\begin{itemize} 1561\item{\tt sic perturbative} This is the default option for the sic 1562directive. After a self-consistent calculation, the Kohn-Sham 1563orbitals are localized with the Foster-Boys algorithm (see section 1564\ref{orbloc}) and the self-interaction energy is added to the total energy. 1565All exchange-correlation functionals implemented in the NWChem can be 1566used with this option. 1567\item{\tt sic oep} With this option the optimized effective potential is 1568built in each step of the self-consistent process. Because the electrostatic 1569potential generated for each orbital involves a numerical 1570integration, this method can be expensive. 1571\item{\tt sic oep-loc} 1572This option is similar to the oep option with the 1573addition of localization of the Kohn-Sham orbitals in each step of the 1574self-consistent process. 1575\end{itemize} 1576With oep and oep-loc options a {\bf xfine grid} (see section \ref{grgrid}) 1577must be used in order to avoid numerical noise, furthermore the hybrid 1578functionals can not be used with these options. More details of the 1579implementation of this method can be found in 1580J. Garza, J. A. Nichols and D. A. Dixon, J. Chem. Phys. 112, 7880 (2000). 1581The components of the sic energy can be printed out using: 1582 1583\begin{verbatim} 1584print "SIC information" 1585\end{verbatim} 1586 1587 1588\section{{\tt MULLIKEN} --- Mulliken analysis} 1589Mulliken analysis of the charge distribution is invoked by the keyword: 1590\begin{verbatim} 1591 MULLIKEN 1592\end{verbatim} 1593When this keyword is encountered, Mulliken analysis of both the input 1594density as well as the output density will occur. 1595For example, to perform a mulliken analysis and print the 1596explicit population analysis of the basis functions, use 1597the following 1598\begin{verbatim} 1599dft 1600 mulliken 1601 print "mulliken ao" 1602end 1603task dft 1604\end{verbatim} 1605 1606 1607\section{{\tt BSSE} --- Basis Set Superposition Error} 1608 1609Particular care is required to compute BSSE by the counter-poise 1610method for the DFT module. In order to include terms deriving from 1611the numerical grid used in the XC integration, the user must label 1612the ghost atoms not just {\tt bq}, but {\tt bq} followed by the given 1613atomic symbol. For example, the first component needed to compute the 1614BSSE for the water dimer, should be written as follows 1615 1616\begin{verbatim} 1617geometry h2o autosym units au 1618 O 0.00000000 0.00000000 0.22143139 1619 H 1.43042868 0.00000000 -0.88572555 1620 H -1.43042868 0.00000000 -0.88572555 1621 bqH 0.71521434 0.00000000 -0.33214708 1622 bqH -0.71521434 0.00000000 -0.33214708 1623 bqO 0.00000000 0.00000000 -0.88572555 1624end 1625 1626basis 1627 H library aug-cc-pvdz 1628 O library aug-cc-pvdz 1629 bqH library H aug-cc-pvdz 1630 bqO library O aug-cc-pvdz 1631end 1632\end{verbatim} 1633 1634Please note that the ``ghost'' oxygen atom has been labeled {\tt bqO}, 1635and not just {\tt bq}. 1636 1637\section{{\tt DISP} --- Empirical Long-range Contribution (vdW)} 1638 1639\begin{verbatim} 1640 DISP 1641 [ vdw <real vdw integer default 2]] 1642 [[s6 <real s6 default depends on XC functional>] \ 1643 [ alpha <real s6 default 20.0d0] \ 1644 [ off ] \ 1645\end{verbatim} 1646 1647When systems with high dependence on van der Waals interactions 1648are computed, the dispersion term may be added empirically 1649through long-range contribution 1650DFT-D, i.e.$ E_{DFT-D}=E_{DFT-KS}+E_{disp}$, 1651%[Q. Wu and W. Yang, \textit{J. Chem. Phys.} \textbf{116} 515 (2002)] 1652where: 1653 1654$$ E_{disp}=-s_6\sum^{N_{atom}-1}_{i=1}\sum^{N_{atom}}_{j=i+1} 1655\frac{C_{6}^{ij}}{R_{ij}^{6}} 1656\left( 1+e^{-\alpha (R_{ij}/R_{vdw}-1)} \right)^{-1}$$ 1657 1658In this equation, the $s_6$ term depends in the functional and basis set used, $C_6^{ij}$ 1659is the dispersion coefficient between pairs of atoms. 1660$R_{vdw}$ and $R_{ij}$ are related with van der Waals atom radii and 1661the nucleus distance respectively. The $\alpha$ value contributes to control the corrections 1662at intermediate distances. 1663 1664There are available three ways to compute $C_6^{ij}$: 1665\begin{enumerate} 1666\item $C_6^{ij}= 1667\frac{2(C_6^{i}C_6^{j})^{2/3}(N_{eff i}N_{eff j})^{1/3}} 1668{C_6^{i}(N_{eff i}^2)^{1/3}+(C_6^{i}N_{eff j}^2)^{1/3}}$ 1669where $N_{eff}$ and $C_6$ are obtained from Q. Wu and W. Yang, \textit{J. Chem. Phys.} 1670\textbf{116} 515 (2002) and U. Zimmerli, M Parrinello and P. Koumoutsakos 1671\textit{J. Chem. Phys.} \textbf{120} 2693 (2004). (Use {\tt vdw 0}) 1672 1673\item $C_6^{ij}=2\frac{C_6^{i}C_6^{j}}{C_6^{i}+C_6^{j}}$. 1674See details in 1675S. Grimme \textit{J. Comp. Chem.} \textbf{25} 1463 (2004). 1676% $R_{vdW}$ and $C_6$ 1677% those were obtained from this paper. 1678(Use {\tt vdw 1}) 1679 1680\item $C_6^{ij}=\sqrt{C_6^{i}C_6^{j}}$ 1681See details in 1682S. Grimme \textit{J. Comp. Chem.} \textbf{27}1787 (2006). 1683% $R_{vdW}$ and $C_6$ 1684% those were obtained from this paper. 1685(Use {\tt vdw 2}) 1686 1687\end{enumerate} 1688 1689Note that in each option there is a certain set of $C_6$ and $R_{vdw}$. 1690 1691For options {\tt vdw 1 } and {\tt vdw 2 }, there are $s_6$ values by default for some functionals 1692and triple-zeta plus double polarization basis set (TZV2P): 1693\begin{itemize} 1694\item {\tt vdw 1}. BLYP 1.40, PBE 0.70 and BP86 1.30. 1695\item {\tt vdw 2}. BLYP 1.20, PBE 0.75, BP86 1.05, B3LYP 1.05, Becke97-D 1.25 and TPSS 1.00. 1696\end{itemize} 1697 1698This capability is also supported for energy gradients and Hessian. 1699Is possible to be deactivated with \verb+OFF+. 1700 1701\section{Print Control} 1702\begin{verbatim} 1703 PRINT||NOPRINT 1704\end{verbatim} 1705 1706The \verb+PRINT||NOPRINT+ options control the level of output in the 1707DFT. Please see some examples using this directive in 1708section \ref{sec:DFTsample}, a sample input file. 1709 Known controllable print options are: 1710 1711\begin{table}[htbp] 1712\begin{center} 1713\begin{tabular}{lcc} 1714 {\bf Name} & {\bf Print Level} & {\bf Description} \\ 1715 ``all vector symmetries'' & high & symmetries of all molecular orbitals \\ 1716 ``alpha partner info'' & high & unpaired alpha orbital analysis \\ 1717 ``common'' & debug & dump of common blocks \\ 1718 ``convergence'' & default & convergence of SCF procedure \\ 1719 ``coulomb fit'' & high & fitting electronic charge density \\ 1720 ``dft timings'' & high & \\ 1721 ``final vectors'' & high & \\ 1722 ``final vector symmetries'' & default & symmetries of final molecular orbitals \\ 1723 ``information'' & low & general information \\ 1724 ``initial vectors'' & high & \\ 1725 ``intermediate energy info'' & high & \\ 1726 ``intermediate evals'' & high & intermediate orbital energies \\ 1727 ``intermediate fock matrix'' & high & \\ 1728 ``intermediate overlap'' & high & overlaps between the alpha and beta sets \\ 1729 ``intermediate S2'' & high & values of S2 \\ 1730 ``intermediate vectors'' & high & intermediate molecular orbitals \\ 1731 ``interm vector symm'' & high & symmetries of intermediate orbitals \\ 1732 ``io info'' & debug & reading from and writing to disk \\ 1733 ``kinetic\_energy'' & high & kinetic energy \\ 1734 ``mulliken ao'' & high & mulliken atomic orbital population \\ 1735 1736 ``multipole'' & default & moments of alpha, beta, and nuclear charge densities \\ 1737 ``parameters'' & default & input parameters \\ 1738 ``quadrature'' & high & numerical quadrature \\ 1739 ``schwarz'' & high & integral screening info \& stats at completion\\ 1740 ``screening parameters'' & high & integral accuracies \\ 1741 ``semi-direct info'' & default & semi direct algorithm \\ 1742\end{tabular} 1743\end{center} 1744\caption{DFT Print Control Specifications} 1745\end{table} 1746 1747 1748 1749 1750