#!/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
#         lim@rani.chem.yale.edu

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(/\.fchk$/, "", 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
	printf("0 COMMENT\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
}