Quantum chemistry computations with Molpro (original course material by Dr. U. W. Schmidt)

To learn how to use the program by yourself or maybe extend your knowledge beyond the
scope of this tutorial one can check out the
links "Quick Start" or "User's Manual" on the following site

To start a calculation we basically need four things:

- the structure of the molecule of interest in (cartesian) coordinates
- the overall charge of our system
- the method that will be used (Hartree-Fock (HF), Density Functional Theory (DFT), Moeller-Plesset perturbation theory (MPx) etc.)
- a basis set.

tar xzvf prakt1.tar.gz

A simple input for the electronic structure calculation of H2 is given below:

basis=6-31G r=1.0 ang geometry={angstrom; h1; h2,h1,r} hf

First, the basis (6-31G) using the so-called Pople notation is specified. G stands for "Gaussian", so a GTO basis with 6 gaussian functions for the inner (non-valence) electrons and 4 (3+1) for
the valence electrons is used. Then, the H-H distance and the so-called internal coordinates (bonds,angles, torsions) are specified.Finally, the key word "hf" stands for "do a Hartree-Fock calculation"

Now we start out first electronic structure calculation by typing

molpro h2.com

During the calculation MOLRPO writes an output file that communicates the most important infos to the user. Please extract the following information:

- number of SCF cycles needed for convergence
- the HF energy (what units are being used?)
- the dipole moment

In the next step we want to perform the visualization of the bonding and antibonding molecular
orbitals (MO) of the H2 molecule. For this we have to redo are electronic structure calculation with a new input file:

basis=6-31G r=1.0 ang geometry={angstrom; h1; h2,h1,r} hf {cube,h2 orbital,all }

and execute molpro again:

molpro h2_orbitals.com

The same calculation as in the previous step will be performed, but additionaly the MOs
and the electron density are written to external files. To view the MOs, we will have to use VMD, a molecular visualization program for displaying, animating, and analyzing large biomolecular systems using 3-D graphics and built-in scripting.
Typing

vmd -f *.cube

will load an existing script that loads the relevant MOs.

To inspect the orbitals, go to Graphics, Representations. Select drawing method CPK. Now you see the two H atoms.

Next, Create Rep (which copies the present representation), and select Drawing method isosurface. With the sliders you can decide the isovalue that is plotted. VMD can show the positive or negative values of the orbital expectation values. To see both, copy the representation Create Rep, and select the negative isovalue points by manualy canging the sign of the Isovalue.

So far, we have only performed so-called single-point calculations, i.e. computed the HF energy
for a fixed nuclear configuration. MOLPRO is quite versatile when it comes to automatically
perform scans along selected coordinates. In our case of the H2, we will now
determine the energy profile (or the potential energy surface (PES)) along the H-H distance. The input
file now looks like this:

r=[0.3 0.4 0.6 0.8 1.0 1.2 1.4 1.6 2.0 4.0 6.0] ang basis=6-31G i=1 geometry={ h1; h2,h1,r(i)} do i=1,#r !loop over range of bond distances hf e1(i)=energy enddo {table,r,e1 title,HF PES for H-H plot,file='h2_hf_pes.plot' }

the program will loop over all specified distances and perform for each of those an independent SCF calculation. Finally, a table of the respective HF energies will be written to a file called *h2_hf_pes.plot*

Now we going to use the plot program "xmgrace" that comes with most Linux distributions
to look at the PES. Typing

xmgrace h2_hf_pes.plot

will reveal the electronic HF energy as a function of the H-H distance.
What is the energy one would expect for the dissociation limit (i.e. r-> oo) of H-H?

As pointed out in the lecture, the HF method gives molecular interaction energies that are only
approximations to the exact quantum many-body energy. In the following we will study the effect
of improving the HF approximation. This will be achieved via perturbation theory, the so-called
Moeller-Plesset (MPx, where x gives the order in the perturbation expansion) and variational
multi-determinant methods, like e.g. configuration interaction (CI). In molpro this can simply be
achieved by adding the respective keywords:

r=[0.3 0.4 0.6 0.8 1.0 1.2 1.4 1.6 2.0 3.0 4.0 5.0 6.0] ang basis=6-31G i=1 geometry={ h1; h2,h1,r(i)} do i=1,#r !loop over range of bond distances hf ;e_hf(i)=energy mp2,NOCHECK ;e_mp2(i)=energy casscf ;e_cas(i)=energy mrci ;e_ci(i)=energy enddo {table,r,e_hf,e_mp2,e_cas,e_ci title,HF PES for H-H plot,file='h2_methods_pes.plot' }

After a sucessful finish of the MOLPRO run you can have a look at the four PESs again using "xmgrace".
In our example the CI method (mrci) can be considered as the most accurate. CASSCF stands for "complete active space SCF", which can be regarded as a truncated CI expansion. Compare your result with respect to the obtained equilibrium bonding distance (How large is the error for this quantity?) and the dissociation limit for
each method. Hereby the CI results should serve as the reference energy curve. Repeat the calculation with changing the
order in the perturbation expansion from 2 to 4 (mp2 -> mp4) and compare the two results. What is the error in the CI result
compared to the analytical results of -1 Hartree. How can we improve the CI dissociation energy?

3 remark line O 0.00000 0.00000 0.00000 H 0.00000 0.00000 0.94000 H 0.00000 0.99030 -0.23192

With this structure, which does not represent the well-known equilibrium structure of water, we first have to perform a so-called geometry optimization. During the geometry optimization cycle a SCF calculation is performed. This one is used to compute the gradient of the total energy on the nuclei, which is then used
to perform a move of the nuclei. Then the new nuclei positions are feeded into a new SCF calculation. These
steps are repeated until the structure (=positions of the nuclei) does not change anymore.

The input file for a geometry optimization looks like this:

basis=6-31G !basis set geomtyp=xyz GEOMETRY=h2o.xyz hf !Hartree-Fock optg,savexyz=h2o_opt.xyz

Compare the optimized geometry "h2o_opt.xyz" with your input geometry "h2o.xyz" using VMD. Determine the bond distances and the bending angle in the water molecule. How does it compare to the experimental results? How can we further improve the structure? How many geometry optimization cycles were necessary to converge the geometry? (you could also plot the number of cycles against the corresponding HF energy using "xmgrace" or "gnuplot")
Using the optimized geometry, we are now computing the MOs of water.
Have a look at the input file "h2o_orbitals.com" first and try to analyze what is going on
in this type of calculation.

basis=6-31G !basis set geomtyp=xyz GEOMETRY=h2o_opt.xyz hf !do scf {cube,h2o orbital,all }

Then type

molpro h2o_orbitals.com

The generated MOs can again be visualized using VMD:

vmd -f *.cube

Have a look at all occupied MOs and elucidate their role for bonding in the water monomer.
Which atomic orbitals of the individual atoms are predominately participating in each MO?
Compare the MOs with the ones shown on the following web site:
Molecular Orbitals of the Water Monomer

Here are the cartesian coordinates:

24 HF-SCF000/STO-3G ENERGY=-667.73565398 C 0.7948047874 -0.2673475369 0.0000000000 C -0.2830459813 -1.0890658393 0.0000000000 N 1.9143454786 -1.1136247133 0.0000000000 N 0.0798949129 -2.4417816138 0.0000000000 N -1.6029576700 -0.6005693362 0.0000000000 C 1.4073520392 -2.3930182886 0.0000000000 C 0.7125900943 1.2091306058 0.0000000000 C -1.8326022679 0.8073510174 0.0000000000 N -0.6606787360 1.6419581534 0.0000000000 O 1.6414281440 2.0047021485 0.0000000000 O -2.9562106712 1.2845856383 0.0000000000 C -0.9238244722 3.0918864981 0.0000000000 C -2.7666488438 -1.4999441112 0.0000000000 C 3.3316990743 -0.7109956220 0.0000000000 H 2.0612180469 -3.2586428364 0.0000000000 H -1.4903555584 3.3842737788 -0.8850934873 H -1.4903555584 3.3842737788 0.8850934873 H 0.0364850129 3.6028877898 0.0000000000 H -3.3823201570 -1.3336563127 -0.8857759769 H -3.3823201570 -1.3336563127 0.8857759769 H -2.4015120487 -2.5249243967 0.0000000000 H 3.9471620995 -1.6091240819 0.0000000000 H 3.5561814227 -0.1173493821 0.8876452774 H 3.5561814227 -0.1173493821 -0.8876452774

and here is what the input file should look like:

basis=6-31G ! basis geomtyp=xyz GEOMETRY=caffeine.xyz set,charge=+0 gdirect hf {cube,caff orbital,homo }

Type to start the HF caclulation:

molpro caffeine.com & tail -f caffeine.out

and observe the information that is written to the screen.
How long does one SCF cycle take on average? How many electrons are involved during this computation?
The so-called HOMO (highest occupied molecular orbital) can be visualized using VMD:

vmd -f caff_orbital_10.2.cube

How many occpuied MO are there?
If there is some time left, try to compute all orbitals by changing "HOMO" to "ALL" in the input file.
Then visualize them along the same lines as in the previous example.