Assignment 1, due date: February 10, 2010 (Dr. Piotr Piecuch) Dipole moment function for the HF molecule (Gaussian and GAMESS calculations). ------------------------------------------------------------------------------ The main purpose of this assignment is the calculation of the dipole moment function for the HF molecule (this function is important in calculations of the infrared absorption line strengths). In general, the dipole moment function describes the dependence of the electronic dipole moment of a molecule on the nuclear geometry. We would like to compare the results of calculations performed with the following three methods: the restricted Hartree-Fock approach (RHF), the coupled-cluster singles and doubles method (CCSD, in two different ways), and the full configuration interaction (FCI) approach. The FCI approach gives the exact dipole moment function for a given basis set. This will be your reference for the analysis of the results. The RHF approach does not describe many-electron correlation effects. These are described by the CCSD method (and FCI, of course). You will be able to learn what the role of the electron correlation effects in describing the dipole moment function of the HF molecule is, and experience the significance of electron correlation in describing bond breaking. ------------------------------------------------------------------------------- Instructions: GAUSSIAN 98 RUNS 1. Prepare inputs for the CCSD (coupled-cluster singles and doubles) calculations for the HF molecule, as described by the DZ basis set. The calculations are for the dipole moment of HF as a function of the internuclear distance R. Here is the set of R values that we want to obtain the dipole moments for: 1.12632, 1.2996, 1.7328 (the approximate equilibrium bond length), 2.166, 2.5992, 3.0324, 3.4656, 4.332, and 5.1984 (all in bohr). The Gaussian input decks should have the following lines: # rccsd=(rw,conver=9)/gen UNITS=Bohr SCF=(QC,Conver=9) # Polar=EnOnly for R < 4.332 bohr, and # rccsd=(rw,conver=8)/gen UNITS=Bohr SCF=(QC,Conver=8) # Polar=EnOnly for R = 4.332 and 5.1984 bohr. In describing the nuclear geometry, use the Cartesian coordinates rather than the Z matrix, orienting the HF molecule along the z axis in the following way: F 0.0 0.0 0.0 H 0.0 0.0 R for example, F 0.0 0.0 0.0 H 0.0 0.0 1.7328 for R = 1.7328 bohr. This is important, since dipole moment is a vector whose components depend on the choice of the molecular frame. The DZ basis set of Dunning which we will use has the following form: H 0 S 3 1.924060D+01 3.282800D-02 2.899200D+00 2.312080D-01 6.534000D-01 8.172380D-01 S 1 1.776000D-01 1.000000D+00 **** F 0 S 6 9.994790D+03 2.017000D-03 1.506030D+03 1.529500D-02 3.502690D+02 7.311000D-02 1.040530D+02 2.464200D-01 3.484320D+01 6.125930D-01 4.368800D+00 2.424890D-01 S 1 1.221640D+01 1.000000D+00 S 1 1.207800D+00 1.000000D+00 S 1 3.634000D-01 1.000000D+00 P 4 4.435550D+01 2.086800D-02 1.008200D+01 1.300920D-01 2.995900D+00 3.962190D-01 9.383000D-01 6.203680D-01 P 1 2.733000D-01 1.000000D+00 **** Make sure that all electrons are correlated in the CCSD runs, which means that you should add the following window at the end of your input file (below the last ****): (...) **** 1,0 An example of an input file for the CCSD calculation for R = 1.7328 bohr will be provided to you during the lab on Wednesday, February 3, 2010 or you can copy it yourself to your directory by typing cp ~piecuch/cem888_SS10/students1/hfccsd1.7328.com ./ 2. Prepare similar input files for the RHF calculations. The Gaussian input deck should begin with the line: # rhf/gen UNITS=Bohr SCF=(QC,Conver=10) An example of an input file for the RHF calculation for R = 1.7328 bohr will be provided to you during the lab on Wednesday, February 3, 2010 or you can copy it yourself to your directory by typing cp ~piecuch/cem888_SS10/students1/hfrhf1.7328.com ./ 3. Once you are done with the calculations, read the values of the dipole moment from the output files. In each case, go to the end of the output file and read the value of the dipole reported as \Dipole= ... You will see three numbers there corresponding to the x, y, and z components of the dipole (z is the internuclear axis). Use the z component for plotting and reporting. The dipole moment values are in atomic units (a.u. or bohr-electron; 1 a.u. = 2.541766 Debye). Use atomic units for reporting and plotting. In these units, the experimental value of the dipole moment (the vibrationless value) is 0.707 a.u. GAMESS RUNS: 1. Prepare inputs for the CCSD property calculations for the HF molecule, as described by the DZ basis set. Again, the calculations are for the dipole moment of HF as a function of the internuclear distance R using the same set of R values as used in the Gaussian calculations. The GAMESS input files should have the following format $CONTRL SCFTYP=RHF CCTYP=CCSD RUNTYP=ENERGY UNITS=BOHR $END $SYSTEM TIMLIM=9999999 MWORDS=20 $END $BASIS GBASIS=DZDUN EXTFIL=.TRUE. $END $CCINP CCPRP=.TRUE. NCORE=0 $END $DATA HF Cnv 2 F 9.0 0.0 0.0 0.0 H 1.0 0.0 0.0 R $END for R = 1.12632, 1.2996, 1.7328, 2.166, 2.5992, 3.0324, 3.4656, and 4.332 bohr. For R = 5.1984 bohr, add the following additional line: $SCF DIIS=.TRUE. SOSCF=.FALSE. $END Note that as in the case of Gaussian calculations, all electrons are correlated (NCORE=0; the default in GAMESS is to freeze core electrons). Use the same DZ basis set as in the Gaussian calculations (called DZDUN in the above input example). The easiest thing to do (the above input format reflects on this) is to place the basis set in a separate file in the same directory as your GAMESS inputs are (in this way, you can build your own library of basis sets, in addition to the basis sets stored in the GAMESS source code, and keep inputs short). In order for GAMESS to read the basis set informaton from the external file (EXTFIL=.TRUE. in the $BASIS input group), you must submit the Gamess runs with the -b flag in gmssub, for example, 'gmsub -b basis hfccsd1.7328', where 'basis' is the name of your file containing the basis set information. You must make sure that the basis set is formatted using the GAMESS style, which is slightly different than the Gaussian style. In the case of the DZ basis set that interests us here, you must use something like this: H DZDUN S 3 1 19.24060000 0.3282800000E-01 2 2.899200000 0.2312080000 3 0.6534000000 0.8172380000 S 1 1 0.1776000000 1.000000000 for the hydrogen atom and F DZDUN S 6 1 9994.790000 0.2017000000E-02 2 1506.030000 0.1529500000E-01 3 350.2690000 0.7311000000E-01 4 104.0530000 0.2464200000 5 34.84320000 0.6125930000 6 4.368800000 0.2424890000 S 1 1 12.21640000 1.000000000 S 1 1 1.207800000 1.000000000 S 1 1 0.3634000000 1.000000000 P 4 1 44.35550000 0.2086800000E-01 2 10.08200000 0.1300920000 3 2.995900000 0.3962190000 4 0.9383000000 0.6203680000 P 1 1 0.2733000000 1.000000000 for fluorine. An example of an input file for the CCSD calculation for R = 1.7328 bohr will be provided to you during the lab on Wednesday, February 3, 2010 or you can copy it yourself to your directory by typing cp ~piecuch/cem888_SS10/students1/hfccsd1.7328.inp ./ 2. Once you are done with the calculations, read the values of the dipole moment from the output files. The RHF dipoles, in Debye, are listed a few lines above --------------------------- COUPLED CLUSTER CALCULATION --------------------------- The CCSD dipoles can be found in two places, a few lines below the line GROUND STATE CCSD PROPERTIES in a.u., and and at the end of the output file in Debye. Use atomic units for reporting and plotting. As mentioned earlier, the unit conversion factor is 1 a.u. = 2.541766 Debye. THINGS TO DO: 1. Report the results by preparing the following table: R RHF(Gaussian) RHF(Gamess) CCSD(Gaussian) CCSD(Gamess) FCI -------------------------------------------- 1.12632 ... ... ... ... ... ... etc. where, as explained earlier, FCI stands for the full configuration interaction approach (which is the exact result for the DZ basis set). The FCI/DZ values of the HF dipole moment (the z component; the x and y components are zero) are: R (bohr) Dipole moment (a.u.) 1.12632 0.75852 1.2996 0.80163 1.7328 0.89789 2.166 0.94752 2.5992 0.92262 3.0324 0.80517 3.4656 0.60494 4.332 0.20937 5.1984 0.05006 2. Prepare a plot with the dipole moment functions obtained with the RHF, CCSD, and FCI methods (use the FCI results provided above). Do it twice, for the GAMESS and Gaussian calculations. You can use gnuplot or any other plotting program you like. A legend and labels for axes are important. Do not forget about them. 3. Comment on the results. Why is the RHF curve so different than FCI or CCSD? What may, in your view, cause the wrong behavior of the RHF dipole moment function, whose values do not approach zero as R approaches infinity (as they should, see the FCI results)? The simple answer "lack of correlation effects in RHF" will not suffice. Try to think about the unphysical components of the RHF wave function that are no longer present when the exact, FCI approach is employed. The dipole moment of HF should vanish at R=infinity, since in this limit the dipole moment of HF becomes a sum of dipole moments of the H and F atoms. 4. Are the dipole moments values obtained in the RHF Gaussian and GAMESS calculations identical? Are the dipole moments values obtained in the CCSD Gaussian and GAMESS calculations identical? Do not worry about the differences at the level of 10^-7 a.u. that are related to convergence criteria used in the calculations, which were not the tightest possible, and focus on the first six decimal places. This is a little tricky, but can you try to explain the differences, if any, based on the contents of your output files and the way the Gaussian and GAMESS dipole moment calculations are executed (outputs are rather transparent about it, but you can always consult the relevant manuals if need be). 5. Examine the Gamess outputs and look for the largest T1 and largest T2 cluster amplitudes printed after the CCSD iterations, which reflect on the role of the leading singly and doubly excited configurations in the wave function expansion relative to the RHF reference determinant. Focus on the largest cluster amplitude in each category and create a table of the values of these amplitudes as a function of the H-F internuclear distance R. Can you determine what orbital excitations are involved in these amplitudes (using the language of sigma, pi, etc. molecular orbitals of the HF molecule)? The outputs have all the information. Can you describe and explain the observed dependence of the largest T1 and T2 amplitudes on the H-F distance R? Please elaborate.