Practical calculation of the B-matrix

To calculate the B-matrix we use the program of Schachtschneider. There are two versions

bmat - no symmetry is included

umat - the point group symmetry of a molecule is included

For the example of water we must include the C2v symmetry of the molecule. We will use the output from the 6-31G HF calculation above to generate the B-matrix. We have written programs to generate the input data file in correct format for bmat and umat. These programs are

bmat_setup - generates an input file for bmat

umat_setup - generates an input file for umat

umat_setup

Before running the program umat_setup you should sketch the molecule and number the atoms in the order that they appear in the Gaussian98 calculation. This numbering must be adhered to in the umat_setup program for consistency. Our sketch for H2O looks like

To run umat_setup we type

Note that the program names are all lower case in the UNIX operating system. The program prompts the user for information. I have included the lines of code that request information and the line FORTRAN code that associates that information with a variable inside the program. The dark red color represents the response that the user will give for this example.

print *, 'Enter the name of the molecule'

read (5,*) molecule

The name of the molecule is a string that will appear in the title of the UMAT input file. Here we input H2O.

print *, 'Enter the point group of the molecule'

read (5,*) pg

The point group of the molecule is a string that will appear in the second line of UMAT input file. This information is used to look up the character table for the point group. Here we input C2V. This tells the UMAT program to use the CN.DAT file for the symmetry information.

print *, 'Enter the Cartesian Coordinate filename'

read (5,*) fname

The Cartesian coordinate file is h2o_631g.coord as generated above from the Standard orientation of the Gaussian98 output file.

print *, 'Enter number of atoms: '

read (5,*) iatoms

3

print *, 'Enter the number of force constants'

read (5,*) nfcs

This variable is used only for the refinement of the F-matrix using the program FFIT by McIntosh and Peterson. At present it is irrelevant and can be set to zero. However, note below:

Since you must work out the group theory for your molecule any way in order to use this program you should have an estimate. For H2O, there are two A1 and one B2 irreducible representations. The A1 representation has two force constants (the symmetric stretch and the bend), and one interaction force constant (total of 3). The B2 representation has 1 force constant. Therefore, nfcs = 4.

print *, 'Enter the number of internal coordinates'

read (5,*) nob

For H2O, this is 3. If there are no redundant coordinates in the problem then the number of force constants will be equal to 3N - 6. If there are redundant coordinates, then this number will be larger. The UMAT program provides you with the option of forming linear combinations of these coordinates to eliminate redundant coordinates. Therefore, the variable nob can be larger than 3N - 6.

Next we enter a loop over the number of internal coordinates. For each internal coordinate we are presented a menu of choices:

print *, 'Enter the number of internal coordinate', i

print *, 'BOND STRETCH - 1'

print *, 'VALENCE ANGLE BEND - 2'

print *, 'OUT-OF-PLANE WAG - 3'

print *, 'TORSION - 4'

print *, 'LINEAR BEND 1 - 5'

print *, 'LINEAR BEND 2 - 6'

read (5,*) code

For a stretch we enter 1.

print *, 'Enter the coordinate label'

read (5,*) coord(i)

The coordinate label is a character string. It contains the atom type and number. Here we enter O1-H2.

type='BOND STRETCH'

print *, 'Enter the atom numbers involved in the ',type

read (5,*) a(i),b(i)

Here we enter

1

2

Note that the order is the same as that in the coordinate label.

print *, 'Enter symmetry factor'

read (5,*) fa(i)

For now just enter 1.

This procedure is repeated for the second stretch and for the valence angle bend. For the bend it is important to enter the atoms in the correct order. From the sketch above we can see that 2 1 3 is one possible order, although 3 1 2 is also valid. The important point is that the pivot or apex atom is atom 1 and this atom is in the middle of the sequence of numbers.

So, for example

print *, 'Enter the coordinate label'

read (5,*) coord(i)

The coordinate label is a character string. It contains the atom type and number. Here we enter H2-O1-H3. Note the numbering scheme above.

type='VALENCE ANGLE BEND'

print *, 'Enter the atom numbers involved in the ',type

read (5,*) a(i),b(i),c(i)

Here we enter

2

1

3

Note that the order is the same as that in the coordinate label.

print *, 'Enter symmetry factor'

read (5,*) fa(i)

Again this is 1.

The torsion and out-of-plane wag coordinates are discussed in a sample problem, ethylene.

Once you have cycled through all of the internal coordinates the umat_setup program will prompt you for an output name.

print *,'Enter output filename: '

read (5,*), fname

h2o_631g.uin

The uin extension means that this is an input file for the program umat.

Now we will run umat.

The first request from umat is:

WHAT IS THE NAME OF THE INPUT DATA FILE?

h2o_631g.uin

WHAT IS THE NAME OF THE OUTPUT FILE?

h2o_631g.uout

WHAT IS THE NAME OF THE SYMMETRY INPUT DATA FILE?

[path]CN.DAT

WHAT IS THE NAME OF THE UMAT OUTPUT FILE?

h2o_631g.umat

The files tmp.ubm and tmp.ufc will be generated automatically when the program is run. The ubm file contains intermediate data for the subroutine UBMAT. The ufc file contains information on the symmetrizing matrices used. In this case the information is obtained from the file CN.DAT.

The file h2o_631g.umat will contain the B-matrix. The output file will have some header lines on top. These should be deleted before this file is used as input for fcartp.

For example, if the file has the appearance

 

WATER ANALYZED USING THE C2V SYMMETRY GROUP

UMAT PROGRAM INPUT

3 9

.0000000000000000E+00 .0000000000000000E+00 .7953159681235821E+00

.0000000000000000E+00 .5846949013903853E+00 -.3976579840617910E+00

.0000000000000000E+00 -.5846949013903853E+00 -.3976579840617910E+00

.0000000000000000E+00 .0000000000000000E+00 -.1741872280742784E+01

.0000000000000000E+00 .5923340685938466E+00 .8709361403713921E+00

.0000000000000000E+00 -.5923340685938466E+00 .8709361403713921E+00

.0000000000000000E+00 -.1169389802780771E+01 .0000000000000000E+00

.0000000000000000E+00 .5846949013903853E+00 -.3976579840617910E+00

.0000000000000000E+00 .5846949013903853E+00 .3976579840617910E+00

 

the lines in yellow should be deleted. In the end the file should contain a line with the number of atoms and number of Cartesian coordinates and then rows of three double precision numbers equal to the product of those numbers.

 

bmat_setup

For larger molecules bmat_setup should be used. An example can be seen in the uracil calculation.