36.24 Exercise: SCF program

Write a closed-shell SCF program for H$_2$O using MATROP!

Hints:

First generate a starting orbital guess by finding the eigenvectors of h0. Store the orbitals in a record. Basis and geometry are defined in the usual way before the first call to MATROP.

Then use a MOLPRO DO loop and call MATROP for each iteration. Save the current energy in a variable (note that the nuclear energy is stored in variable ENUC). Also, compute the dipole moment in each iteration. At the end of the iteration perform a convergence test on the energy change using the IF command. This must be done outside MATROP just before the ENDDO. At this stage, you can also store the iteration numbers, energies, and dipole moments in arrays, and print these after reaching convergence using TABLE. For the following geometry and basis set

geometry={o;h1,o,r;h2,o,r,h1,theta}   !Z-matrix geometry input
r=1 ang                               !bond length
theta=104                             !bond angle
basis=vdz                             !basis set
thresh=1.d-8                          !convergence threshold

the result could look as follows:

 SCF has converged in 24 iterations 

    ITER       E             DIP
     1.0   -68.92227207    2.17407361
     2.0   -71.31376891   -5.06209922
     3.0   -73.73536433    2.10199751
     4.0   -74.64753557   -1.79658706
     5.0   -75.41652680    1.43669203
     6.0   -75.77903293    0.17616098
     7.0   -75.93094231    1.05644998
     8.0   -75.98812258    0.63401784
     9.0   -76.00939154    0.91637513
    10.0   -76.01708679    0.76319435
    11.0   -76.01988143    0.86107911
    12.0   -76.02088864    0.80513445
    13.0   -76.02125263    0.83990621
    14.0   -76.02138387    0.81956198
    15.0   -76.02143124    0.83202128
    16.0   -76.02144833    0.82464809
    17.0   -76.02145450    0.82912805
    18.0   -76.02145672    0.82646089
    19.0   -76.02145752    0.82807428
    20.0   -76.02145781    0.82711046
    21.0   -76.02145792    0.82769196
    22.0   -76.02145796    0.82734386
    23.0   -76.02145797    0.82755355
    24.0   -76.02145797    0.82742787

It does not converge terribly fast, but it works!