1% $Id$ 2\label{sec:interface} 3 4NWChem has interfaces to several different packages which are listed below. 5In general, the NWChem authors work with the authors of the other packages 6to make sure that the interface works. However, any problems with the 7interface should be reported to the 8{\tt nwchem-users@emsl.pnl.gov} e-mail list. 9 10 11\section{{\tt DIRDYVTST} --- DIRect Dynamics for Variational Transition State Theory} 12\label{sec:dirdyvtst} 13 14by Bruce C. Garrett,\\ 15Environmental Molecular Sciences Laboratory,\\ 16Pacific Northwest Laboratory, Richland, Washington\\ 17\\ 18Yao-Yuan Chuang and Donald G. Truhlar,\\ 19Department of Chemistry and Super Computer Institute,\\ 20University of Minnesota, MN 55455-0431\\ 21\\ 22and interfaced to NWChem by\\ 23\\ 24Ricky A. Kendall,\\ 25Scalable Computing Laboratory,\\ 26Ames Laboratory and Iowa State University, Ames, IA 50011\\ 27\\ 28Theresa L. Windus,\\ 29Environmental Molecular Sciences Laboratory,\\ 30Pacific Northwest Laboratory, Richland, Washington 31 32If you use the DIRDYVTST portion of NWChem, please use following citation 33in addition to the usual NWChem citation from Section \ref{sec:intro}: 34\begin{quote} 35 DIRDYVTST, Yao-Yuan Chuang and Donald G. Truhlar, 36 Department of Chemistry and Super Computer Institute, 37 University of Minnesota; Ricky A. Kendall,Scalable Computing Laboratory, 38 Ames Laboratory and Iowa State University; Bruce C. Garrett and Theresa L. 39 Windus, Environmental Molecular Sciences Laboratory, Pacific Northwest 40 Laboratory. 41\end{quote} 42 43\subsection{Introduction} 44 45By using DIRDYVTST, a user can carry out electronic structure calculations 46with NWChem and use the resulting energies, gradients, and Hessians for 47direct dynamics calculations with \htmladdnormallink{POLYRATE} 48{http://comp.chem.umn.edu/polyrate/}. 49This program prepares the file30 input for \htmladdnormallink{POLYRATE} 50{http://comp.chem.umn.edu/polyrate/} from NWChem 51electronic structure calculations of energies, gradients and Hessians at the 52reactant, product, and saddle point geometries and along the minimum 53energy path. Cartesian geometries for the reactants, products, and 54saddle points need to be input to this program; optimization of 55geometries is not performed in this program. Note that \verb+DIRDYVTST+ is 56based on the \htmladdnormallink{DIRDYGAUSS program} 57{http://comp.chem.umn.edu/dirdyg/} 58and is similar to two other programs: \htmladdnormallink{DDUTILITIES} 59{http://comp.chem.umn.edu/ddutil/} and \htmladdnormallink{GAUSSRATE} 60{http://comp.chem.umn.edu/gaussrate/}. Users of this module are 61encouraged to read the 62\htmladdnormallink{POLYRATE} 63{http://comp.chem.umn.edu/polyrate/} manual since they will need to 64create the file fu5 input to run calculations with \htmladdnormallink{POLYRATE} 65{http://comp.chem.umn.edu/polyrate/}. 66 67Notes about the code: 68 69Input. The code has been written to parallel, as much as possible, 70the \verb+POLYRATE+ code. 71 72Output. There is one default output file for each \verb+DIRDYVTST+ run - .file30. 73 74Integrators for following the reaction path. 75Currently the Euler and three Page-McIver (PM) methods 76are implemented. The PM methods are the local quadratic approximation 77(LQA), the corrected LQA (CLQA), and the cubic (CUBE) algorithm. 78The PM methods are implemented so that the Hessian can be reused at 79intermediate steps at which only the gradient is updated. 80 81 82\subsection{Files} 83 84Test runs are located in directories in \verb+$NWCHEM_TOP/QA/tests+. Test 85runs are available for two systems: $H + H_2$ and $OH + H_2$. 86 87The $H + H_2$ test uses the Euler integration method at the SCF/3-21G level 88of theory to calculate points along the reaction path. 89This test is located in the \verb+$NWCHEM_TOP/QA/tests/h3tr1+ directory. 90 91The $OH + H_2$ test uses the Page-McIver CUBE algorithm to calculate points 92on the SCF/3-21G surface and does additional single point calculations at 93the SCF/6-31G* level of theory. This test is located in the 94\verb+$NWCHEM_TOP/QA/tests/oh3tr3+ directory. 95 96Note: These tests are set up with SCF, however, other levels of 97theory can be used. The initial hessian calculations at the reactants, 98products and saddle point can cause some problems when numerical hessians 99are required (especially when there is symmetry breaking in the wavefunction). 100 101\subsection{Detailed description of the input} 102 103The input consists of keywords for NWChem and keywords related to \verb+POLYRATE+ 104input. The first set of inputs are for NWChem with the general input block 105of the form: 106 107\begin{verbatim} 108 DIRDYVTST [autosym [real tol default 1d-2] | noautosym] 109 [THEORY <string theory> [basis <string basis default "ao basis">] \ 110 [ecp <string ecp>] [input <string input>]] 111 [SPTHEORY <string theory> [basis <string basis default "ao basis">] \ 112 [ecp <string ecp>] [input <string input>]] 113 ... 114 END 115\end{verbatim} 116 117\subsubsection{Use of symmetry} 118The use of symmetry in the calculation is controlled by the keyword 119\verb+autosym | noautosym+ which is used as described in the geometry 120directive (see Section \ref{sec:geom}). 121\verb+Autosym+ is on by default. A couple words of warning here. 122The tolerance related to \verb+autosym+ can cause problems when taking the 123initial step off of the transition state. If the tolerance is too large and 124the initial step relatively small, 125the resulting geometry will be close to a higher symmetry than is really 126wanted and the molecule will be symmetrized into the higher symmetry. 127To check this, the code prints out the symmetry at each geometry along the 128path. It is up to the user to check the symmetry and make sure that 129it is the required one. In preverse cases, the user may need to turn 130\verb+autosym+ off (\verb+noautosym+) if changing the tolerance doesn't 131produce the desired results. In the case that autosym is used, the 132user does not need to worry about the different alignment of the molecule 133between NWChem and POLYRATE, this is taken care of internally in the 134DIRDYVTST module. 135 136\subsubsection{Basis specification} 137The basis name on the theory or sptheory directive is that 138specified on a basis set directive (see Section \ref{sec:basis}) and 139{\em not} the name of a standard basis in the library. If not 140specified, the basis set for the sptheory defaults to the 141theory basis which defaults to \verb+"ao basis"+. 142 143\subsubsection{Effective core potentials} 144If an effective core potential is specified in the usual fashion (see 145Section \ref{sec:ecp}) outside of the \verb+DIRDYVTST+ input then this will be 146used in all calculations. If an alternative ECP name (the name 147specified on the ECP directive in the same manner as done for basis 148sets) is specified on one of the theory directives, then this ECP will 149be used in preference for that level of theory. 150 151\subsubsection{General input strings} 152 153For many purposes, the ability to specify the theory, basis and 154effective core potential is adequate. All of the options for each 155theory are determined from their independent input blocks. However, 156if the same theory (e.g., DFT) is to be used with different options 157for theory and sptheory, then the general input strings must 158be used. These strings are processed as NWChem input each time the 159theoretical calculation is invoked. The strings may contain any NWChem 160input, except for options pertaining to \verb+DIRDYVTST+ and the task directive. 161The intent is that the strings be used just to control the options 162pertaining to the theory being used. 163 164A word of caution. Be sure to check that the options are producing 165the desired results. Since the NWChem database is persistent, 166the input strings should fully define the calculation you wish to have happen. 167 168For instance, if the theory model is DFT/LDA/3-21g and the 169sptheory model is DFT/B3LYP/6-311g**, the \verb+DIRDYVTST+ input might look like this 170\begin{verbatim} 171 dirdyvtst 172 theory dft basis 3-21g input "dft\; xc\; end" 173 sptheory dft basis 6-311g** input "dft\; xc b3lyp\; end" 174 .... 175 end 176\end{verbatim} 177The empty \verb+XC+ directive restores the default LDA 178exchange-correlation option (see Section \ref{sec:xc}). Note that 179semi-colons and other quotation marks inside the input string must be 180preceded by a backslash to avoid special interpretation. 181 182\subsubsection{POLYRATE related options} 183 184These keyword options are simlar to the \verb+POLYRATE+ input format, except there 185are no ENERGETICS, OPTIMIZATION, SECOND, TUNNELING, and RATE sections. 186 187\paragraph{GENERAL section} 188 189 190The GENERAL section has the following format: 191 192\begin{verbatim} 193*GENERAL 194 [TITLE <string title>] 195 ATOMS 196 <integer num> <string tag> [<real mass>] 197 ... 198 END 199 [SINGLEPOINT] 200 [SAVEFILE (vecs || hess || spc) 201\end{verbatim} 202 203 Descriptions 204 205 TITLE is a keyword that allows the user to input a description of 206 the calculation. In this version, the user can only have a 207 single-line description. 208 209 For example: 210 TITLE Calculation of D + HCl reaction 211 212 ATOMS is a list keyword that is used to input a list of 213 the atoms. It is similar to \verb+POLYRATE+ in that the order of the 214 atom and the atomic symbol are required in a single line. If 215 isotope of the element is considered then the atomic mass is 216 required in units of amu. 217 218 For example: 219\begin{verbatim} 220 ATOMS 221 1 H 2.014 222 2 H 223 3 Cl 224 END 225\end{verbatim} 226 227 SINGLEPOINT is a keyword that specifies that a single 228 point calculation is to be performed at the reactants, 229 products and saddle point geometries. The type of 230 single point calculation is specified in the \verb+sptheory+ 231 line. 232 233 SAVEFILE is a keyword that specifies that NWChem files 234 are to be saved. Allowed values of variable input to 235 SAVEFILE are vecs, hess, and spc for saving the files 236 base theory movecs, base theory hessian and singlepoint 237 calculation movecs. 238 239\paragraph{REACT1, REACT2, PROD1, PROD2, and START sections} 240 241These sections have the following format: 242\begin{verbatim} 243*(REACT1 || REACT2 || PROD1 || PROD2 || START) 244 GEOM 245 <integer num> <real x y z> 246 ... 247 END 248 SPECIES (ATOMIC || LINRP || NONLINRP || LINTS || NONLINTS default NONLINRP) 249\end{verbatim} 250 251 REACT1 and REACT2 are input for each of the reactants and PROD1 and PROD2 252 are input for each of the products. REACT1 and PROD1 are required. START 253 is the input for the transition state if one exists, or starting point to 254 follow downhill the MEP. 255 256 Descriptions 257 258 GEOM is a list keyword that indicates the geometry of the molecule 259 in Cartesian coordinates with atomic unit. 260 261 For example: 262\begin{verbatim} 263 GEOM 264 1 0.0 0.0 0.0 265 2 0.0 0.0 1.5 266 END 267\end{verbatim} 268 269 SPECIES is a variable keyword that indicates the type of the 270 molecule. Options are: ATOMIC (atomic reactant or product), 271 LINRP (linear reactant or product), NONLINRP 272 (nonlinear reactant or product), LINTS (linear transition 273 state), and NONLINTS (nonlinear transition state). 274 275 For example: 276 SPECIES atomic 277 278\paragraph{PATH section} 279 280The Path section has the format: 281 282\begin{verbatim} 283*PATH 284 [SCALEMASS <real scalemass default 1.0>] 285 [SSTEP <real sstep default 0.01>] 286 [SSAVE <real ssave default 0.1>] 287 [SHESS <real shess default SSAVE>] 288 [SLP <real slp default 1.0>] 289 [SLM <real slm default -1.0>] 290 [SIGN (REACTANT || PRODUCT default REACTANT)] 291 [INTEGRA (EULER || LQA || CLQA || CUBE default EULER)] 292 [PRINTFREQ (on || off default off)] 293\end{verbatim} 294 295 Descriptions 296 297 SCALEMASS is a variable keyword that indicates the arbitrary 298 mass (in amu) used for mass-scaled Cartesian coordinates. 299 This is the variable called mu in published papers. Normally, 300 this is taken as either 1.0 amu or, for bimolecular reactions, 301 as the reduced mass of relative translation of the reactants. 302 303 SSTEP is a variable keyword that indicates the numerical step 304 size (in bohrs) for the gradient grid. This is the step 305 size for following the minimum energy path. 306 307 SSAVE is a variable keyword that indicates the numerical step 308 size (in bohrs) for saving the Hessian grid. At each 309 save point the potential and its first and second 310 derivatives are recalculated and written to the .file30 311 file. For example, if SSTEP=0.01 and SSAVE=0.1, then the 312 potential information is written to .file30 every 10 313 steps along the gradient grid. 314 315 SHESS is a variable keyword that indicates the numerical step 316 size (in bohrs) for recomputing the Hessian when using a 317 Page-McIver integrator (e.g., LQA, CLQA, or CUBE). For 318 Euler integration SHESS = SSAVE. For intermediate points 319 along the gradient grid, the Hessian matrix from the 320 last Hessian calculation is reused. For example, if 321 SSTEP=0.01 and SHESS=0.05, then the Hessian matrix is 322 recomputed every 5 steps along the gradient grid. 323 324 SLP is a variable keyword that indicates the positive limit of 325 the reaction coordinate (in bohrs). 326 327 SLM is a variable keyword that indicates the negative limit of 328 the reaction coordinate (in bohrs). 329 330 SIGN is a variable keyword used to ensure the conventional definition 331 of the sign of s, $s < 0$ for the reactant side and 332 $s > 0$ for the product side, is followed. \verb+PRODUCT+ 333 should be used if the eigenvector at the saddle point points 334 toward the product side and \verb+REACTANT+ if the 335 eigenvector points toward the reactant side. 336 337 INTEGRA is a variable keyword that indicates the integration 338 method used to follow the reaction path. Options are: EULER, 339 LQA, CLQA, and CUBE. 340 341 PRINTFREQ is a variable keyword that indicates that projected 342 frequencies and eigenvectors will be printed along the MEP. 343 344\subsubsection{Restart} 345 346\verb+DIRDYVTST+ calculations should be restarted through the normal NWChem 347mechanism (See Section \ref{sec:start}). The user needs to change the 348\verb+start+ directive to a \verb+restart+ directive and get rid of any 349information that will overwrite important information in the RTDB. The 350\verb+file.db+ and \verb+file.file30+ need to be available for the 351calculation to restart properly. 352 353\subsubsection{Example} 354 355This is an example that creates the file30 file for POLYRATE for $H + H_2$. 356Note that the multiplicity is that of the entire supermolecule, a doublet. 357In this example, the initial energies, gradients, and Hessians are calculated 358at the UHF/3-21G level of theory and the singlepoint calculations are 359calculated at the MP2/cc-pVDZ level of theory with a tighter convergence 360threshold than the first SCF. 361 362\begin{verbatim} 363start h3test 364 365basis 366 h library 3-21G 367end 368 369basis singlepoint 370 h library cc-pVDZ 371end 372 373scf 374 uhf 375 doublet 376 thresh 1.0e-6 377end 378 379dirdyvtst autosym 0.001 380 theory scf input "scf\; uhf\; doublet\; thresh 1.0e-06\; end" 381 sptheory mp2 basis singlepoint input \ 382 "scf\; uhf\; doublet\; thresh 1.0e-07\; end" 383*GENERAL 384 TITLE 385 Test run: H+H2 reaction, Page-McIver CLQA algorithm, no restart 386 387 ATOMS 388 1 H 389 2 H 390 3 H 391 END 392 393 SINGLEPOINT 394 395*REACT1 396 GEOM 397 1 0.0 0.0 0.0 398 2 0.0 0.0 1.3886144 399 END 400 401 SPECIES LINRP 402 403*REACT2 404 GEOM 405 3 0.0 0.0 190.3612132 406 END 407 408 SPECIES ATOMIC 409 410*PROD2 411 GEOM 412 1 0.0 0.0 190.3612132 413 END 414 415 SPECIES ATOMIC 416 417*PROD1 418 419 GEOM 420 2 0.0 0.0 1.3886144 421 3 0.0 0.0 0.0 422 END 423 424 SPECIES LINRP 425 426*START 427 428 GEOM 429 1 0.0 0.0 -1.76531973 430 2 0.0 0.0 0.0 431 3 0.0 0.0 1.76531973 432 END 433 434 SPECIES LINTS 435 436*PATH 437 SSTEP 0.05 438 SSAVE 0.05 439 SLP 0.50 440 SLM -0.50 441 SCALEMASS 0.6718993 442 INTEGRA CLQA 443end 444 445task dirdyvtst 446\end{verbatim} 447 448