Quantum chemistry computations with Molpro (original course material by Dr. U. W. Schmidt)
In the previous lecture we've learned about the theoretical basis of electronic
structure methods, and in particular about the Hartree-Fock method.
Now it is time to jump into the electronic structure business and perform
a couple of calculations by yourself.
The package we are going to use is called MOLPRO, which has been developed
by Peter Knowles and Hans-Juergen Werner over the last 15 years. It is among
the modern electronic structure codes available and provides an easy-to-use scripting language that allows for
a user-friendly access to quantum chemistry. Additionaly MOLPRO provides for the largest set of implemented post-Hartree-Fock methods.
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:
tar xzvf prakt1.tar.gz
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:
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
Finally we are gonna to compute the electronic structure of the well-known
Caffeine molecule, which is a xanthine alkaloid compound that acts as a stimulant in humans.
The structure of this biomolecule is shown to the left.
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.