Subject: Creating input files for PSI88 package from Gaussian Formcheck files

From: Dongchul Lim <[email-masked]>
Date: Tue, 26 Jul 94 22:44:36 EDT

***************************************************
***** Update on the PSI88 Package of Programs *****
***************************************************

____________________________________________________________________________

This post is regarding how to generate input files for PSI88 package
>from a Gaussian-related files.
I'm posting this message on the request of one of co-authors of PSI88.
So this is official almost :-).

D. Lim, Department of Chemistry, Yale University.
[email-masked]
____________________________________________________________________________

1. What is 'PSI88' and where can I get it?

        The PSI88 is a MO plotting package, which is written by Prof.
    William Jorgensen and Dr. Severance at Yale University and available from
    www.ccl.net via anonymous ftp. This package is designed to plot 2D
    representation of molecular orbitals from given molcular coordinates and
    orbital coefficients. It supports semi-empirical, STO-3G, 3-21(++)G(**) 
    and 6-31(++)G** basis sets. It consists of three subprograms 'PSI1',
    'PSICON' and 'PSI2'. 'PSI1' is responsible for generation of 3D grid of
    wave function psi (hence it was named 'PSI'). PSICON reads in the 3D grid
    and generates contour lines (isovalued lines in x,y,z directions).
    Finally 'PSI2' reads in these contour lines, performs hidden-line removal
    and gives plots for HPGL, PostScript, etc.

2. What are the typical methods of getting molecular orbital coefficients?

        Any semi-empirical and ab initio package should output orbital
    coefficients. I'm focusing only to the Gaussian ab initio package.
    In calculations with Gaussian, orbital coefficients can be obtained 
    mainly in three ways.

 1) %chk=filename
        If '%chk=filename' is specified in the zmatrix, a checkpoint file
    is created in the calculation. It'll contain all necessay information 
    to generate input files for the PSI88. A utility called 'chk2psi' is
    provided in PSI88 distribution for this purpose.

  2) POP=FULL
	If the keyword 'POP=FULL' is used in the zmatrix, MO cofficients
    are printed in the Gaussian output file. Usually the last occurrenct
    of "Molecular Orbital Coefficients' are the coefficients of the final
    optimized geometry.

  3) FORMCHECK=ALL
	If the keyword 'FORMCHECK=ALL' is used in the zmatrix, a formcheck
    file is created in the calculation. It contains coordinates, orbital
    energies, orbital coefficients, etc. Since this file contains most
    reliable information on molecular orbitals, I recommend anyone to use
    formcheck files for creating PSI88 input files. An nawk script 'fchk2psi'
    (which I'm going to enclose here) can be used for this purpose.

3. What utilities are available for generating input files for the PSI88
   package from Gaussian-related files?

  1) chk2psi
	This is written to generate PSI88 input files from a Gaussian check-
    point file. However, cautions should be taken before its use, i.e.,
    it always extracts the Z-matrix orientation and may give wrong MO plots.
    Note that this comes with the PSI88 distribution.

  2) g2psi
    I wrote a utilty which extracts molecular coordinates and oribtal
    coefficients from a Gaussian output which has been executed with POP=FULL
    option. Interested users may want to contact me to get it mailed to you.

  3) fchk2psi
	This script is written in AWK programming language which is an inter-
    preted language and available in most UNIX systems as /usr/bin/nawk.
    It extracts molecular coordinates and MO coefficients from a formcheck
    file and creates input files for the PSI88 package.
    It will be found at the end of this posting.

    Users of the PSI88 should realize that the 'Standard' orientation of
    molecule must be used for plotting MO's IF IT EXISTS. Otherwise, the
    'Z-matrix orientation' can be used instead. A formcheck file contains
    always the right molecular orientation for MO's. So I recommend everyone
    to use formcheck files for creating PSI88 input files.
    The name of the script is recommended to be 'fchk2psi'.
    The Usage is simply "fchk2psi filename". If filename has a suffix ".log",
    e.g., "file.log", PSI input file names will be "file.psi1", "file.psi2"
    and "file.psicon". Otherwise, "psi1", etc will be simply appended to
    the filename. If input for fchk2psi comes from standard input (stdin),
    e.g., "zcat filename.Z | fchk2psi" or "fchk2psi < filename",
    files named "tmp.psi1", "tmp.psi2" and "tmp.psicon" will be created.
    
____________________________________________________________________________

Using your favorite editor, carve out the script between 'CUT HERE' pairs.
Name it 'fchk2psi' and make it executable by doing 'chmod +x fchk2psi'.
Finally place it in a directory which is included in your evironment
variable 'PATH'.

--------------------- CUT HERE --------------------------------------------
#!/usr/bin/nawk -f
# @(#) fchk2psi: extracts molecular coordinates and MO coefficients 
#      from Gaussian formchekpoint file and creates input files for
#      the PSI88 package written by Prof. William Jorgensen and Dr.
#      Dan Severance at Yale University.
#      To get MO coefficients from Gaussian, specify FORMCHECK=ALL.
#
# Usage:  fchk2psi filename
# Author: Dongchul Lim
#         Department of Chemistry, Yale University
#         [email-masked]

BEGIN	{
	natoms = 0;		# number of atoms
	nbasis = 0;		# number of basis functions
	bohr2A = 0.529177;	# Bohr to Angstrom conversion factor
	number = "^[-+]?([0-9]+[.]?[0-9]*|[.][0-9]+)([eE][-+]?[0-9]+)?$";

	getline;		# 1st line (ignore)
	getline;		# 2nd line (contains basis set)
	basis = substr($0,35,20);	# basis set
	sub(/^[ \t]+/, "", basis);	# remove preceding spaces and tabs

	if (FILENAME == "-") {	# standard input
		file_psi1   = "tmp.psi1";
		file_psi2   = "tmp.psi2";
		file_psicon = "tmp.psicon";
	} else {
		file = FILENAME;
		sub(/\.log$/, "", file);	# replace ".log" by ""
		file_psi1   = sprintf("%s.psi1", file);
		file_psi2   = sprintf("%s.psi2", file);
		file_psicon = sprintf("%s.psicon", file);
	}

if (Debug) {
	printf("basis  = %s\n", basis);
	printf("psi1   = %s\n", file_psi1);
	printf("psi2   = %s\n", file_psi2);
	printf("psicon = %s\n", file_psicon);
}

}

/Number of basis functions/	{
	nbasis = substr($0,50,20) + 0;
}
/Number of atoms/	{
	natoms = substr($0,50,20) + 0;
}
/Atomic numbers/	{
	n = 0;
	while (getline > 0 && match($1,number)) {
		for(i=1;i<=NF;i++) {n++; an[n] = $i+0;}
	}
}
/Current cartesian coordinates/	{
	n = 0;
	while (getline > 0 && match($1,number)) {
		for(i=1;i<=NF;i++) {n++; xyz[n] = $i * bohr2A;}
	}
}
/Alpha MO coefficients/	{
	print_psi1(xyz, an, natoms, basis, file_psi1);
	n = 0;
	while (getline > 0 && match($1,number)) {
		for(i=1;i<=NF;i++) {
			n++;
			mo[n] = $i+0;
			if (n == 8) {
				print_mo(mo,n, file_psi1);
				n = 0;
			}
		}
	}
	print_mo(mo,n, file_psi1);
}

END	{
	if (natoms == 0) {
		printf("\007No atoms found.\n");
		exit;
	}

	# psi1 has been already printed. So output just psi2 and psicon.
	print_psi2(xyz, an, natoms, basis, file_psi2);
	print_psicon(basis, file_psicon);
}

function	print_mo(mo, n,	file,	i)
{
	for(i=1;i<=n;i++) printf("% 10.6f", mo[i]) > file
	if (n>0) printf("\n") > file
}

function	print_psi1 (xyz, an, natoms, basis, file)
{
	printf("%s\n", basis) > file
	printf("AUTO0\n") > file
	printf("0101  1.0\n") > file
	printf("0\n") > file
	print_coord(xyz, an, natoms, file);
	printf("99\n") > file
}

function	print_coord (xyz, an, natoms, file,	i)
{
	for(i=0;i<natoms;i++) {
		printf("%2d        % 10.6f% 10.6f% 10.6f\n",
 			an[i+1], xyz[i*3+1], xyz[i*3+2], xyz[i*3+3]) > file
	}
}

function	print_psicon (basis, file)
{
	printf("%s\n", basis) > file
	printf(" 1 1 0 1\n")  > file
	printf(" 0.075000\n") > file
}

function	print_psi2 (xyz, an, natoms, basis, file)
{
	printf("%s\n", basis)  > file
	printf("SUBTITLE\n")   > file
	printf("010000 1.0\n") > file
	printf("00\n")         > file
	print_coord(xyz, an, natoms, file);
	printf("99\n") > file
	printf("  61.1000  132.1000    1.1000    0.8500\n") > file
	printf("02\n") > file
}

--------------------- CUT HERE --------------------------------------------