First, it should be emphasized that the developers of CFOUR encourage you to use analytic second derivatives whenever possible, especially when one is planning to calculate an anharmonic force field. The reason is simply that the cubic
and quartic force constants that are needed for VPT2 energy levels and intensities (or other applications) are computed
by numerical differentiation. If you base the calculation on analytic second derivatives, this (obviously) means that your
higher-order force constants will be more accurate than if the calculation is based on analytic first derivatives or (please don't do this unless it is absolutely necessary) energies.
With that caveat, it should also be mentioned that a lot of methods (EOM-CC in particular) do not yet have analytic second derivative implementations, and anharmonic force fields via EOM is something that several people want/need to do. Such efforts must - for now - be driven by VIB=FINDIF (analytic first derivatives used to generate second derivatives) calculations. One could create a ZMAT file like:
water
H
O 1 R
H 2 R 1 A
R = 0.964203501252364
A = 103.602079370900668
*CFOUR(CALC=CCSD[T],BASIS=ANO0,CC_CONV=9,LINEQ_CONV=9
SCF_CONV=9,DROPMO=1
ANHARM=VPT2
VIB=FINDIF)
and submit the job to a queue, but the following major consideration arises (pretend that this is a BIG job, like EOM-CCSD with 350 basis functions):
The job might not finish, due to computer crash, lack of convergence in the SCF, lack of convergence in CC equations, power failure, etc.
Thus, it is not unwise to break this job into a number of smaller jobs (one gradient per job), and this can now be done quite easily with scripts provided here. The procedure to carry out a harmonic frequency calculation, using finite differences of gradients with each job run separately is dealt with exhaustively in
Calculating harmonic frequencies by finite differences in parallel
and
Scripts for calculating harmonic frequencies by finite differences in parallel
which you are encouraged to read (and practice with) before running a VPT2 calculation in the manner described here.
With all of the above as an introduction, here is what one should do:
0. Add the following keywords to the ZMAT file: ANH_ALGORITHM=PARALLEL,FREQ_ALGORITHM=PARALLEL and
(this is very important) FD_PROJECT=OFF. That is, the file:
water
H
O 1 R
H 2 R 1 A
R = 0.964203501252364
A = 103.602079370900668
*CFOUR(CALC=CCSD[T],BASIS=ANO0,CC_CONV=9,LINEQ_CONV=9
SCF_CONV=9,DROPMO=1
ANHARM=VPT2
ANH_ALGORITHM=PARALLEL,FREQ_ALGORITHM=PARALLEL,FD_PROJECT=OFF
VIB=FINDIF)
serves as an example.
1. Follow the procedure described in Scripts for calculating harmonic frequencies by finite differences in parallel
2. In the directory where the harmonic calculation is processed (see step #4 in Calculating harmonic frequencies by finite differences in parallel), run the following script:
Script pfindif_proc_1
#!/bin/csh
xjoda
xsymcor
foreach i (fja.*)
cp $i FJOBARC
xja2fja
xsymcor
end
xjoda
rm zmat*
xcubic
which is identical to the script used in step 4 of Calculating harmonic frequencies by finite differences in parallel except that the first set of zmatxxx files are removed and the xcubic executable is now invoked at the end. (You could, of course, strictly adhere to the directions of step #1 above and then just remove the zmat* files and invoke xcubic after you have finished).
3. At this point, you have a directory with a JOBARC and JAINDX file that contain the information from the harmonic job and the zmatxxx jobs for each displaced point used to compute the cubic and quartic force fields. These latter zmatxxx files are similar in spirit to that at the top of this directory in that they are input files for a SECOND derivative calculation
(albeit one by finite difference) and therefore each serves as the template for an individual harmonic frequency calculation at each displaced point. With the JOBARC and JAINDX files in place as well as the various zmatxxx files and the run.dummy script discussed elsewhere (see Scripts for calculating harmonic frequencies by finite differences in parallel), run the following script:
script pfindif_anharm_generate
#!/bin/csh
#!
#! Run this script after doing the harmonic calculation. Should (well, must)
#! be run from the directory where the zmatxxxx files are located, each of which
#! is the file that specifies the harmonic calculation for a displaced point.
#!
foreach i (zmat*)
set j = `echo $i | sed 's/zmat//' `
#! Loop through files and make separate directories for each of the zmat files
#! present
mkdir $j
cp $i $j/ZMAT
ln -s $basisdir/GENBAS $j/GENBAS
#! If you are using ECPDATA uncomment the following line
#! ln -s $ecpdir/ECPDATA $j/ECPDATA
cd $j
#! Go into each directory and generate the individual gradient zmat files and
#! then make runscripts from the master run.dummy file in main directory
xjoda > xjodaout
xsymcor > xsymcorout
foreach k (zmat*)
set l = `echo $k | sed 's/zmat//' `
sed s/job.xxxxx/job.$j.$l/g < ../run.dummy > run.tmp
sed /'set workdir'/d < run.tmp > run.tmp2
sed '/set tmpdir/ i\set workdir=\('$PWD''\) < run.tmp2 > run.tmp3
sed '/set tmpdir/s/xxxxx/'$j.$l'/' < run.tmp3 > run.tmp4
sed s/xxxxx/$l/g < run.tmp4 > run.$j.$l
rm run.t*
end
cd ..
else
endif
end
(note that $basisdir needs to be defined. If you are using ECPDATA then you
have to define $ecpdir also). Using the ZMAT file above, the appearance of the working directory after invoking this script is:
001 DCT FCM000 fja.004 JMOLplot sub3 zmat005
002 DIPDER FCMFINAL fja.005 JOBARC VMLSYM
003 DIPOL FCMINT FJOBARC MOL ZMAT
004 dnor FILES FRQARC MOLDEN zmat001
005 dnorraw fja.001 GENBAS NORMCO zmat002
ALLDONE dvr.in fja.002 GRD QUADRATURE zmat003
coriolis.joda FCM fja.003 JAINDX run.dummy zmat004
Each of the directories 001, 002 etc. contain the input files and scripts needed to run the individual gradients for the displaced point harmonic frequency calculations. For example, here are the contents of 002:
[stanton@hbar 2]$ ls 002
FILES JMOLplot MOLDEN run.002.003 xjodaout zmat001 zmat004
GENBAS JOBARC run.002.001 run.002.004 xsymcorout zmat002 zmat005
JAINDX MOL run.002.002 run.002.005 ZMAT zmat003
4. Run all of these displaced point gradient calculations. Depending on the size of your molecule, there can be lots and lots of these (there are forty-two for our example). All of the files run.xxx.yyy, where ... is the displaced point number for the frequency calculation and yyy is the displaced point number for the gradient calculation, must be run. When they finish, you'll have lots of fja files. Here are the contents of the 002 directory after all of its jobs have finished:
[stanton@hbar 2]$ ls 002
FILES fja.004 JMOLplot out.001 out.005 run.002.004 ZMAT zmat004
fja.001 fja.005 JOBARC out.002 run.002.001 run.002.005 zmat001 zmat005
fja.002 GENBAS MOL out.003 run.002.002 xjodaout zmat002
fja.003 JAINDX MOLDEN out.004 run.002.003 xsymcorout zmat003
5. Everything is in place now, and just requires processing. Return to the main directory (the one in which the harmonic frequencies were processed, and which contains the 001, 002 etc. subdirectories, and run the following script:
script pfindif_anharm_procall
#!/bin/csh
mkdir proc
rm proc\*
echo "Processing initial harmonic frequency calculation"
xclean
xjoda > anharm.out
xsymcor >> anharm.out
foreach i (fja*)
cp $i FJOBARC
xja2fja >> anharm.out
xsymcor >> anharm.out
end
xjoda >> anharm.out
xcubic >> anharm.out
cp JOBARC jobarc.0
cp JAINDX jaindx.0
foreach i (0*)
echo "----------------------------------------------------------------------" >> anharm.out
echo "Processing output from displaceed point $i" >> anharm.out
echo "Processing output from displaceed point $i"
echo "----------------------------------------------------------------------" >> anharm.out
cd $i
xclean
xjoda >> ../anharm.out
xsymcor >> ../anharm.out
foreach j (fja*)
cp $j FJOBARC
xja2fja >> ../anharm.out
xsymcor >> ../anharm.out
end
xjoda >> ../anharm.out
rm FJOBARC
xja2fja >> ../anharm.out
cp FJOBARC ../proc/fja_all.$i
cd ..
end
echo "----------------------------------------------------------------------" >> anharm.out
echo "Processing final output" >> anharm.out
echo "Processing final output"
echo "----------------------------------------------------------------------" >> anharm.out
cd proc
cp ../jobarc.0 JOBARC
cp ../jaindx.0 JAINDX
foreach i (fja*)
cp $i FJOBARC
xja2fja >> ../anharm.out
xcubic >> ../anharm.out
end
6. Inspect the anharm.out file for your output. Hopefully everything has gone well.
Our test example produces the following anharmonic frequencies:
Harmonic Fundamental Anharmonic Harmonic Fundamental Anharm
Mode Frequency Frequency Contribution Intensity Intensity Contrib
7 1659.6193 1610.5745 -49.0447 73.0281 72.7522 -0.2759
8 3834.8937 3653.0864 -181.8073 6.3744 3.0884 -3.2860
9 3958.1028 3760.5259 -197.5769 42.6482 33.3371 -9.3111
Once you have done this a few times, you can just think of it this way:
1. Create your ZMAT file and place a GENBAS file in the directory.
2. Run the pfindif_setup script.
3. Send the jobs that it generates off to the queues.
4. Run the pfindif_proc_1 script.
5. Run the pfindif_anharm_generate script.
6. Run the (very many) jobs that this generates, all of which are in the 0** subdirectories.
7. After they are done, run pfindif_anharm_procall.
And "that's it"!
CFOUR is partially supported by the U.S. National Science Foundation.