File: running.tex

package info (click to toggle)
psicode 3.3.0-3
  • links: PTS, VCS
  • area: main
  • in suites: lenny
  • size: 32,284 kB
  • ctags: 12,083
  • sloc: ansic: 218,247; cpp: 35,679; fortran: 10,489; perl: 5,413; sh: 3,881; makefile: 1,874; yacc: 110; lex: 52; csh: 12
file content (671 lines) | stat: -rw-r--r-- 27,217 bytes parent folder | download | duplicates (2)
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.