1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390 391 392 393 394 395 396 397 398 399 400 401 402 403 404 405 406 407 408 409 410 411 412 413 414 415 416 417 418 419 420 421 422 423 424 425 426 427 428 429 430 431 432 433 434 435 436 437 438 439 440 441 442 443 444 445 446 447 448 449 450 451 452 453 454 455 456 457 458 459 460 461 462 463 464 465 466 467 468 469 470 471 472 473 474 475 476 477 478 479 480 481 482 483 484 485 486 487 488 489 490 491 492 493 494 495 496 497 498 499 500 501 502 503 504 505 506 507 508 509 510 511 512 513 514 515 516 517 518 519 520 521 522 523 524 525 526 527 528 529 530 531 532 533 534 535 536 537 538 539 540 541 542 543 544 545 546 547 548 549 550 551 552 553 554 555 556 557 558 559 560 561 562 563 564 565 566 567 568 569 570 571 572 573 574 575 576 577 578 579 580 581 582 583 584 585 586 587 588 589 590 591 592 593 594 595 596 597 598 599 600 601 602 603 604 605 606 607 608 609 610 611 612 613 614 615 616 617 618 619 620 621 622 623 624 625 626 627 628 629 630 631 632 633 634 635 636 637 638 639 640 641 642 643 644 645 646 647 648 649 650 651 652 653 654 655 656 657 658 659 660 661 662 663 664 665 666 667 668 669 670 671
|
\subsection{Geometry Specification}
Full molecular geometry has to be specified in form of Cartesian coordinates or
a Z-matrix. Cartesian coordinates of atoms are specified via a keyword
\keyword{geometry} which has to be a member of either \keyword{default}
or \keyword{input} sections:
\begin{verbatim}
geometry = (
(atomname1 x1 y1 z1)
(atomname2 x2 y2 z2)
(atomname3 x3 y3 z3)
..........
(atomnameN xN yN zN)
)
\end{verbatim}
where \keyword{atomname$i$} can take the following values:
\begin{itemize}
\item element symbol (H, He, Li, Be, B, etc.);
\item full element name (hydrogen, helium, lithium, etc.);
\item {\em ghost} atom symbol (G) or name (ghost). Ghost atom is an atom
of formal charge 0.0, it can be useful to specify the location of
the off-nucleus basis functions;
\item {\em dummy} atom symbol (X). Dummy atoms can be useful only to specify
Z-matrices of proper symmetry (not used in \PSIthree; see below) or
which contain linear fragments.
\end{itemize}
Hence the following two examples are equivalent to one another:
\begin{verbatim}
geometry = (
(H 0.0 0.0 0.0)
(f 1.0 0.0 0.0)
(Li 3.0 0.0 0.0)
(BE 6.0 0.0 0.0)
)
\end{verbatim}
\begin{verbatim}
geometry = (
(hydrogen 0.0 0.0 0.0)
(FLUORINE 1.0 0.0 0.0)
(Lithium 3.0 0.0 0.0)
(berillium 6.0 0.0 0.0)
)
\end{verbatim}
The keyword \keyword{units} specifies the units for the coordinates:
\begin{itemize}
\item \keyword{units = bohr} -- atomic units (Bohr), default;
\item \keyword{units = angstrom} -- angstroms ($\AA$);
\end{itemize}
The \keyword{zmat} keyword can be used to specify a Z-matrix for the molecule.
It also has to be put in either \keyword{default} or \keyword{input} sections.
The format of this vector is as follows:
\begin{verbatim}
zmat = (
(atomname1)
(atomname2 ref21 bond_dist2)
(atomname3 ref31 bond_dist3 ref32 bond_angle3)
(atomname4 ref41 bond_dist4 ref42 bond_angle4 ref43 tors_angle4)
(atomname5 ref51 bond_dist5 ref52 bond_angle5 ref53 tors_angle5)
...........................
(atomnameN refN1 bond_distN refN2 bond_angleN refN3 tors_angleN)
)
\end{verbatim}
where
\begin{itemize}
\item \keyword{bond\_dist$i$} is the distance (in units specified by
keyword \keyword{units}) from nucleus number $i$ to
nucleus number \keyword{ref$i$1}. The units
\item \keyword{bond\_angle$i$} is the angle formed by nuclei $i$,
\keyword{ref$i$1}, and \keyword{ref$i$2};
\item \keyword{tors\_angle$i$} is the torsion angle formed by nuclei $i$,
\keyword{ref$i$1}, \keyword{ref$i$2}, and \keyword{ref$i$3};
\end{itemize}
Some care has to be taken when constructing a Z-matrix for a molecule
which contains linear fragments. For example, let's construct a Z-matrix
for a linear conformation of HNCO. The first three atoms (HNC) can be specified
as is, but the fourth atom (O) poses a problem -- the torsional angle cannot be
defined with respect to the linear HNC fragment. The solution is to add
2 dummy atoms to the definition:
\begin{verbatim}
zmat = (
(h)
(n 1 1.012)
(x 2 1.000 1 90.0)
(c 2 1.234 3 90.0 1 180.0)
(x 4 1.000 2 90.0 3 180.0)
(o 4 1.114 5 90.0 2 180.0)
)
\end{verbatim}
Of course, a choice of the method for geometry specification is solely
a matter of convenience.
\subsection{Molecular Symmetry}
\PSIthree\ can determine automatically the largest Abelian point group
for a valid framework of centers (the framework also includes ghost
atoms, but it does not include dummy atoms).
It will then use the symmetry porperties of the system in computing energy,
forces, and other properties.
However, in certain instances it is desirable to use a lower than the
full symmetry. The keyword \keyword{subgroup} is used to specify a subgroup of
the full molecular point group. The allowed values are \keyword{c2v},
\keyword{c2h}, \keyword{d2}, \keyword{c2}, \keyword{cs}, \keyword{ci},
and \keyword{c1}. For certain combinations of a group and
its subgroup there is no unique way to determine which subgroup is
implied. For example, $D_{\rm 2h}$ has 3 non-equivalent $C_{\rm 2v}$ subgroups,
e.g. $C_{\rm 2v}(X)$ consists of symmetry operations $\hat{E}$, $\hat{C}_2(x)$,
$\hat{\sigma}_{xy}$, and $\hat{\sigma_{xz}}$.
To specify subgroups precisely one has to use the \keyword{unique\_axis}
keyword. E.g. the following input will specify the $C_{\rm 2v}(X)$
subgroup of $D_{\rm 2h}$ to be the computational point group:
\begin{verbatim}
input: (
geometry = (
........
)
units = angstrom
subgroup = c2v
unique_axis = x
)
\end{verbatim}
\subsection{Basis Sets}
An atomic basis set is normally identified by a
string. Currently, there exist three ways to specify which basis sets
to use for which atoms:
\begin{itemize}
\item \keyword{basis = string} -- all atoms use basis set type. If the
basis string contains any ``special'' characters (e.g., parentheses,
asterisks) then the string must be enclosed in quotation marks, e.g.,
"6-311++G(d,p)".
\item \keyword{basis = (string1 string2 string3 ... stringN)} --
\keyword{string {\em i}}
specifies the basis set for atom {\em i}. Thus, the number of strings
in the \keyword{basis} vector has to be the same as the number of
atoms (including ghost atoms but excluding dummy atoms). Another
restriction is that symmetry equivalent atoms should have same basis
sets, otherwise \PSIinput\ will use the string provided for the
so-called unique atom out of the set of symmetry equivalent ones.
\item
\begin{verbatim}
basis = (
(element1 string1)
(element2 string2)
...........
(elementN stringN)
)
\end{verbatim}
\keyword{string {\em i}} specifies the basis set for chemical element
\keyword{element {\em i}}.
\end{itemize}
\subsubsection{Default Basis Sets}
\PSIthree\ default basis sets are located in \pbasisdat\ which is located in
{\tt \$psipath/share}. Table \ref{table:basisset} lists basis sets
defined in \pbasisdat.
\begin{table}[tbp]
%\special{landscape}
\caption{~~~Basis sets available in PSI 3.2}
%\vspace{0.015in}
\begin{center}
\begin{tabular}{|l|l|l|}
\hline
\hline
\hline
Basis Set &Atoms &Aliases\\
\hline
\hline
\textbf{Huzinaga-Dunning} & &\\
\hline
(4S/2S) & H &\\
(9S5P/4S2P) & B-F &\\
(11S7P/6S4P) & Al-Cl &\\
DZ & H, B-F, Al-Cl &\\
DZP & H, B-F, Al-Cl &\\
DZ-DIF & H, B-F, Al-Cl &\\
DZP-DIF & H, B-F, Al-Cl &\\
\hline
\hline
\textbf{Wachters} & &\\
\hline
WACHTERS & K, Sc-Cu &\\
WACHTERS-F & Sc-Cu &\\
\hline
\hline
\textbf{Pople-type} & &\\
\hline
STO-3G & H-Ar &\\
3-21G & H-Ar &\\
6-31G & H-Ar, K, Ca, Cu &\\
6-31G* & H-Ar, K, Ca, Cu &6-31G(d)\\
6-31+G* & H-Ar &6-31+G(d)\\
6-31G** & H-Ar, K, Ca, Cu &6-31G(d,p)\\
6-311G & H-Ar &\\
6-311G* & H-Ar &6-311G(d)\\
6-311+G* & H-Ne &6-311+G(d)\\
6-311G** & H-Ar &6-311G(d,p)\\
6-311G(2df,2pd) & H-Ne &\\
6-311++G** & H, B-Ar &6-311++G(d,p)\\
6-311G(2d,2p) & H-Ar &\\
6-311++G(2d,2p) & H-Ar &\\
6-311++G(3df,3pd) & H-Ar &\\
\hline
\hline
\textbf{Triple-Zeta} & &\\
\hline
TZ2P & H, B-F, Al-Cl &\\
TZ2PD & H &\\
TZ2PF & H, B-F, Al-Cl &\\
TZ-DIF & H, B-F, Al-Cl &\\
TZ2P-DIF & H, B-F, Al-Cl &\\
TZ2PD-DIF & H &\\
TZ2PF-DIF & H, B-F, Al-Cl &\\
\hline
\hline
\textbf{Correlation Consistent} & &\\
\textbf{ (N = D,T,Q,5,6)} & & \\
\hline
CC-PVNZ & H-Ar &cc-pVNZ\\
CC-PV(N+D)Z & Al-Ar &cc-pV(N+d)Z\\
CC-PCVNZ & B-Ne &cc-pCVNZ\\
AUG-CC-PVNZ & H-He, B-Ne, Al-Ar &aug-cc-pVNZ\\
AUG-CC-PV(N+D)Z & Al-Ar &aug-cc-pCV(N+d)Z\\
AUG-CC-PCVNZ (N${<}$6) & B-F &aug-cc-pCVNZ\\
D-AUG-CC-PVNZ & H &\\
PV7Z & H, C, N, O, F, S &pV7Z\\
AUG-PV7Z & H, C, N, O, F, S &aug-pV7Z\\
AUG-CC-PV7Z & H, N, O, F &aug-cc-pV7Z\\
\hline
\hline
\hline
\end{tabular}
\end{center}
\end{table}
\subsubsection{Custom Basis Sets}
To make a custom basis set, enter the information in either of the
following four files:
\begin{itemize}
\item \pbasisdat\ -- only if you think it should be added to \PSIthree.
You also might want to check in your additions and changes so that
everyone could benefit from them. Refer to the \PSIthree\ Programmer's
Manual for information on how to access \PSIthree\ repository.
\item An arbitrary text file. To specify the file's location use keyword
\keyword{basisfile}:
\begin{verbatim}
input: (
% The meaning of this is pretty obvious
basisfile = "/home/users/tool/chem/h2o/mybasis.in"
% If the location ends with '/', "basis.dat" is automatically appended
% Hence this specifies /home/users/tool/chem/basis.dat !
basisfile = "/home/users/tool/chem/"
)
\end{verbatim}
Use this option if you want to use the basis set file in a project
which involves running more than one computation.
\item File named \basisdat, which resides in the working directory along
with \inputdat . Same use as the previous entry.
\item \inputdat
\end{itemize}
The order in which \PSIinput\ program searches for basis sets is the order
in which files appear in our checklist.
A contracted Cartesian Gaussian-Type Orbital
\begin{eqnarray}
\phi_{\rm CGTO} & = & x^ly^mz^n\sum_i C_i \exp(-\alpha_i[x^2+y^2+z^2])
\end{eqnarray}
where
\begin{eqnarray}
L & = & l+m+n
\end{eqnarray}
is written as
\begin{verbatim}
basis: (
ATOM_NAME: "BASIS_SET_LABEL" = (
(L (C1 alpha1)
(C2 alpha2)
(C3 alpha3)
...
(Cn alpha4))
)
)
\end{verbatim}
To scale a basis set, a scale factor may be added as the last item
in the specification of each contracted Gaussian function. For
example, to scale the S functions in a 6-31G** basis for hydrogen,
one would use the following
\begin{verbatim}
hydrogen:"6-31G**" =
( (S ( 18.73113696 0.03349460)
( 2.82539437 0.23472695)
( 0.64012169 0.81375733) 1.2 )
(S ( 0.16127776 1.00000000) 1.2 )
(P ( 1.10000000 1.00000000))
)
\end{verbatim}
In this example, both contracted S functions have their exponents
scaled by a factor of (1.2)$^2$ = 1.44. The output file should show
the exponents after scaling.
\subsection{Electronic Structure Specification}
The reference electronic configuration of a molecule is specified via
a combination of keywords \keyword{reference} and
\keyword{multiplicity} and occupation vectors \keyword{docc} and
\keyword{socc}. However, the latter may not be necessary as \PSIcscf\
may guess occupations (\keyword{docc} and \keyword{socc} arrays) for
you had \keyword{charge}, \keyword{multiplicity}, and
\keyword{reference} have been specified. It is the easiest way to
specify electronic configuration for your system, but remember that
guessing algorithms \PSIcscf\ uses are far from perfect. Hence you
should check guessed occupations every time you let \PSIcscf\ guess
them for you.
To determine the electronic occupations in \PSIthree\ manually, first
construct symmetry orbitals using group theory and fill them according
to regular valence bond arguments. To define your occupations in
\PSIthree, use the \keyword{docc} and \keyword{socc} arrays. But only
\keyword{docc} and \keyword{socc} may not be enough to specify
precisely the spin couplings in your system. That's where
\keyword{reference} and \keyword{multiplicity} keywords come
in. \keyword{multiplicity} is equal 2S+1, where S is the spin quantum
number of the system. \keyword{reference} can equal
\begin{itemize}
\item rhf (default) - spin-restricted reference for closed shell molecules.
\keyword{multiplicity} may only equal to 1 in this case.
\item rohf - spin-restricted reference for open shell molecules.
If multiplicity=1 and socc has two singly occupied orbitals in
different symmetry blocks - it's equivalent to the old
opentype=singlet statement. Otherwise it's assumed to be a
high-spin open-shell case (equivalent to the old
opentype=highspin statement).
\item uhf - spin-unrestricted reference for closed shell or
high-spin (parallel spins) open shell system.
\item twocon for two determinantal wavefunctions. The largest
component should be specified by the docc and socc arrays.
Multiplicity has to be set to 1.
\end{itemize}
For $^1{\rm A}_1$ methylene, the occupation is
(1a1)2(2a1)2(1b2)1(3a1)2 so the docc is:
\begin{verbatim}
reference = rhf or uhf
multiplicity = 1
docc = (3 0 0 1)
\end{verbatim}
For the $^3{\rm B}_1$ state of methylene, the electronic configuration
is (1a1)2(2a1)2(1b2)2(3a1)1(1b1)1 so the docc and socc arrays are:
\begin{verbatim}
reference = rohf or uhf
multiplicity = 3
docc = (2 0 0 1)
socc = (1 0 1 0)
\end{verbatim}
For the $^1{\rm B}_1$ state of methylene however, the docc and socc
arrays are also:
\begin{verbatim}
reference = rohf
multiplicity = 1
docc = (2 0 0 1)
socc = (1 0 1 0)
\end{verbatim}
Since most of the basis sets are highly contracted in the core regions,
core electrons are routinely frozen and corresponding virtual
orbitals are deleted. This is accomplished via the \keyword{frozen\_docc}
and \keyword{frozen\_uocc} arrays. Simply specify the symmetry of
the frozen orbital and \PSIthree\ will do the rest.
To freeze the lowest $a_1$ orbital and delete the corresponding
highest $b_2$ orbital, they would look like this:
\begin{verbatim}
frozen_docc = (1 0 0 0)
frozen_uocc = (0 0 0 1)
\end{verbatim}
\subsection{Single-Point Energy Computation}
Along with the wavefunction type, the nuclear framework, basis set,
and electronic configuration are sufficient
for a single-point evaluation of the electronic energy.
The electronic wavefunction is specified via the \keyword{wfn} keyword.
The range of allowed wavefunctions is listed in Table \ref{table:methods}.
\subsection{Geometry Optimization}
\PSIoptking\ is the program responsible for orchestrating the process
of geometry optimization. It can do a number of tasks automatically,
such as generating internal coordinates, produce empirical force
constant matrix, if necessary, update it, and check if geometry
optimization is over. Some or all of the following files are necessary
to perform a geometry optimization with \PSIoptking:
\begin{itemize}
\item \FILE{11.dat} - contains the cartesian geometry and the nuclear
forces, produced by \PSIcderiv ;
\item \fconstdat - contains force
constants; if absent - empirical force constants will be generated by
\PSIoptking ;
\item \intcodat - contains internal coordinates in a
format readable by a human; if absent - internal coordinates are
generated automatically by \PSIoptking .
\end{itemize}
The procedure for setting up such a calculation is as follows:
\begin{itemize}
\item define internal coordinates if desired (or, do nothing, and
\PSIoptking\ will do it for you automatically!)
\item obtain a set of force constants in an fconst.dat file (or,
again, \PSIoptking\ can do this automatically for you)
\item If analytic gradients are available for your chosen method, set
\keyword{dertype=first}. If not, set \keyword{dertype=none} and
also set \keyword{numerical\_dertype=first} in the
\keyword{default} section.
\item Run the optimization by setting the \keyword{opt} flag set to true
and \keyword{ nopt} to the number of geometry optimization
steps (say, around 5 to 10). If analytic gradients are not
available, then \keyword{nopt} instead gives the number of energy
points to compute. This should be the desired number of
geometry optimization steps, multiplied by (2*num\_symm\_coord + 1),
where num\_symm\_coord is the number of totally-symmetric internal
coordinates.
\end{itemize}
The precision with which geometry is optimized depends on the residual
forces on the nuclei. By default \PSIoptking\ will terminate the job
if the residual cartesian gradients in \FILE{11.dat} are less than
1E-5 in atomic units. It is probably enough for most
tasks. Going below this will most likely waste CPU
time unless you are doing benchmarks.
An important aspect of a geometry optimization is the accuracy of the
first derivatives of energy that \PSIthree\ computes. Depending on
how poorly your wavefunction has been convereged, the gradients
themselves may not be sufficiently accurate for the requested
convergence criterion. After computing first derivatives of the
energy, \PSIcints\ runs a simple check of the quality of the energy
derivative. It's a good idea to look at \PSIcints ' output to make
sure that the gradients are OK.
Let us take a look at each step involved in optimizing molecular geometry.
\subsubsection{Internal Coordinates and Structure of \keyword{intco} Vector}
This section is largely obsolete now with the addition of the \PSIoptking\
program which can generate internal coordinates automatically. At present
\PSIoptking\ cannot handle molecules larger than a few
atoms but it should change in the immediate future. Hence
you may still specify internal coordinates manually
as described here, but this ability may become obsolete someday.
\PSIthree\ currently carries out all optimizations in internal
coordinates. The internals are specified in either \inputdat\ or
\intcodat. First, the primitive internals are defined. These are
individual stretches, bends, torsion, out-of-plane deformations, and
two different linear bends denoted lin1 and lin2. All of these are
defined in Wilson, Decius, and Cross. An example for methane is below:
\begin{verbatim}
intco: (
stre = (
(1 1 2)
(2 1 3)
(3 1 4)
(4 1 5)
)
bend = (
(5 2 1 5)
(6 3 1 5)
(7 4 1 5)
(8 2 1 4)
(9 3 1 4)
(10 2 1 3)
)
\end{verbatim}
After the primitives are defined, they are constructed into
symmetrized internals with the totally symmetric placed in
the SYMM vector and the rest placed in the ASYMM vector. For
optimizations, only the SYMM internals need to be defined.
Likewise, if during an optimization a molecule breaks symmetry,
the internals have been improperly defined. Again, methane is done
below:
\begin{verbatim}
symm = (
("(1) stretch"(1 2 3 4))
)
asymm = (
("(2) E bend"(10 7))
("(3) T2 stretch"(1 2 -3 -4 ))
("(4) T2 bend"(10 -7))
("(5) E torsion"(8 -5 -9 6))
("(6) T2 bend" (8 -6))
("(7) T2 bend"( 5 -9))
("(8) T2 stretch"(1 3 -2 -4))
("(9) T2 stretch"(1 4 -2 -3))
)
)
\end{verbatim}
The SYMM and ASYMM vectors have two or three components: the first
is a label enclosed by quotation marks and the second is the list
of primitive internals comprising this vector. Some
internals have been multiplied by -1 to reflect the appropriate
symmetries. If the internals need to be weighted by some prefactor,
then a third vector may be used:
\begin{verbatim}
symm = (
("generic coord" (1 -2 -3) (2.0 1.0 1.0))
)
\end{verbatim}
For more information in defining symmetric internals, refer to Cotton's text.
\subsubsection{Force Constant Matrix and Structure of \fconstdat}
The quality of the force constants, or Hessian, is critical for
optimizing weakly bound structures. In order to start an optimization,
one needs the \fconstdat\ file. For those of you that can speak
Fortran 77, this file is written in 8F10.7 format. It is the lower
triangle of the force constant matrix in internal coordinates. The
order of the forces is identical to the order of the SYMM and ASYMM
vectors. For the methane-water dimer, an excerpt from a real
\fconstdat\ is shown below: \begin{verbatim}
5.654908
.217027 5.616085
-.006096 -.001145 6.078154
-.055291 .026921 .023732 .317485
.004063 -.155489 -.003880 -.170201 .873999
-.146900 -.285108 .001153 -.023346 .196008 .605239
.001037 -.001959 .002945 -.024597 .014432 .005302 .388622
\end{verbatim}
Ideally, your diagonal elements should be the much larger than the
non-diagonal elements. If you need an \fconstdat\ file, you have four
options:
\begin{enumerate}
\item Create a diagonal matrix of 1's
\item Create a diagonal matrix with 5 for stretching coordinates,
2 for bending coordinates, and 1 for all other coordinates
\item Let \PSIoptking\ generate an empirical Hessian for you
\item Run a second derivative to obtain a Cartesian Hessian
and transform that to internals( fconst.dat) with intder.
\end{enumerate}
Clearly, the list starts at the most
approximate and gets more accurate.
\subsection{Frequency Analysis}
Currently, analytic second derivatives are avaliable only for closed-shell
Hartree-Fock (the preliminary implementation is not yet very efficient).
For methods with analytic gradients, one must carry out finite difference
calculations to obtain second derivatives. If one has the
misfortune of needing frequencies for methods where only energies
exist(e.g. FCI or MP2-R12), finite differences with gradients
approximated by energy points must be used. At present, obtaining the
Hessian numerically from energies or gradients is not automated by the
program. There are two types of frequency analysis programs available
within \PSIthree, \PSInormco\ and \PSIintder. For more details
regarding these programs, see the manual pages for each program
respectively.
\subsection{Evaluation of one-electron properties}
For now, take a look at the available documentation for \PSIoeprop .
\subsection{Plotting one-electron properties}
Program \PSIoeprop\ can evaluate certain one-electron properties on a
grid of points and then print them out in format suitable for image
rendering with external programs. Currently, electron and spin
density, electron and spin density gradient, Laplacian of electron and
spin density, and electrostatic potential can be evaluated on an
arbitrary rectangular two-dimensional grid and output in a format
suitable for feeding to a interactive visualization program {\tt
PlotMTV} (version 1.3 and higher). {\tt PlotMTV} is a freeware code
developed by Kenny Toh. It can be downloaded off many web sites in
source or binary form.
In addition, values of molecular orbitals can be evaluated on an
arbitrary rectangular three-dimensional grid and output for furher
rendering of high-quality images with a program {\tt MegaPov} (version
0.5). {\tt MegaPov} is an unofficial patch for a ray-tracing code
{\tt POV-Ray}. Information on {\tt MegaPov} can be found at
\htmladdnormallink{{http://nathan.kopp.com/patched.htm}}{http://nathan.kopp.com/patched.htm}.
Let's look at how to set up input for spin density evaluation on a
two-dimensional grid. An input secion of \PSIoeprop\ might look like
this:
\begin{verbatim}
oeprop:(
grid = 2
spin_prop = true
grid_origin = (0.0 -5.0 -5.0)
grid_unit_x = (0.0 1.0 0.0)
grid_unit_y = (0.0 0.0 1.0)
grid_xy0 = (0.0 0.0)
grid_xy1 = (10.0 10.0)
nix = 20
niy = 20
)
\end{verbatim}
\keyword{grid} specifies the type of a property and the type of a grid
\PSIoeprop\ needs to compute. Allowed values are 0 (default, no
property evaluation on a grid), 1 (electrostatic potential on a 2D
grid), 2 (electron/spin density on a 2D grid), 3 (gradient of
electron/spin density on a 2D grid), 4 (Laplacian of electron/spin
density on a 2D grid), 5 (molecular orbital values on a 3D
grid). Since \keyword{spin\_prop}\ is set, the spin density will be
evaluated on a grid.
Grid specification is a little bit tricky but very
flexible. \keyword{grid\_origin}\ specifies the origin of the
rectangular coordinate system associated with the grid in the
refernence frame. \keyword{grid\_unit\_x}\ specifies a reference frame
vector which designates the direction of the x-axis of the grid
coordinate system. \keyword{grid\_unit\_y}\ is analogously a
reference frame vector which, along with the \keyword{grid\_unit\_x},
completely specifies the grid coordinate system.
\keyword{grid\_unit\_x}\ and \keyword{grid\_unit\_y}\ do not have to
be normalized, neither they need to be orthogonal to either other -
orthogonalization is done automatically to ensure that unit vectors of
the grid coordinate system are normalized in the reference frame too.
\keyword{grid\_xy0}\ is a vector in the grid coordinate system that
specifies a vertex of the grid rectangle with the most negative
coordinates. Similarly, \keyword{grid\_xy1}\ specifies a vertex of the
the grid rectangle diagonally opposite to \keyword{grid\_xy0}.
Finally, \keyword{nix}\ and \keyword{niy}\ specify the number of
intervals into which the $x$ and $y$ sides of the grid rectangle are
subdivided. To summarize, the above input specifies a rectangular (in
fact, square) 21 by 21 grid of dimensions 10.0 by 10.0 lying in the
$yz$ plane and centered at origin.
Running \PSIoeprop\ on such input will create a file called
\file{sdens.dat} (for file names refer to man page on \PSIoeprop),
which can be fed directly to {\tt PlotMTV} to plot the data.
Specification of a three-dimensional grid for plotting MO isosurfaces
({\tt grid = 5}) is just slightly more complicated. The index of the
MO which needs to be plotted is specified by keyword
\keyword{mo\_to\_plot}. The index is specified in Pitzer order and
not according to the orbital energies. The reference frame is
specified by keywords \keyword{grid\_origin}, \keyword{grid\_unit\_x}\
and \keyword{grid\_origin\_y}\ (the third axis of the grid coordinate
system is specified by by the vector product of
\keyword{grid\_unit\_x}\ and \keyword{grid\_unit\_y}). Since in this
case we are dealing with the three-dimensional grid coordinate system,
one needs to specify two diagonally opposite vertices of the grid box
via \keyword{grid\_xyz0}\ and \keyword{grid\_xyz1}. The number of
intervals along $z$ is specified via \keyword{niz}. The final input
may look like this:
\begin{verbatim}
oeprop:(
grid = 5
mo_to_plot = 10
grid_origin = (-5.0 -5.0 -5.0)
grid_unit_x = (1.0 0.0 0.0)
grid_unit_y = (0.0 1.0 0.0)
grid_xyz0 = (0.0 0.0 0.0)
grid_xyz1 = (10.0 10.0 10.0)
nix = 20
niy = 20
niz = 20
)
\end{verbatim}
Running \PSIoeprop\ on input like this will produce two files
\begin{itemize} \item \file{mo.dat} - MO values tabulated in a column
in order where $x$ runs fast; \item \file{mo.pov} - the {\tt MegaPov}
command file which sets a number of parameters. You might have to
adjust a number of parameters manually to obtain the best looking
output. \end{itemize} Rendering of an image can be done via {\tt
megapovplus +Imo.pov}. Consult {\tt MegaPov} documentation for more
information on how to use it.
\subsection{Utilities}
\subsubsection{\PSIgeom}
The program \PSIgeom\ reads a set of Cartesian coordinates and
determines from them the bond distances (Bohr and angstrom), bond
angles, torsional angles, out-of-plane angles (optional), moments of
inertia, and rotational constants. It requires either a \FILE{11.dat}
or \geomdat\ and writes \geomout.
|