#!/bin/sh
# @(#) g2psi: extracts coordinates and MO coefficients 
# and generates psi88 input files (i.e., psi1, psi2 and psicon).
# Molecular coordinates are extracted from the last occurrence
# of "Standard orientation". If the "Standard orientation" does
# not exist, "Z-Matrix orientation" is extracted instead.
# Molecule orbital coefficients are extracted from the last
# occurrence of "Molecular Orbital Coefficients".

printUsage()	{
	echo "Usage: `basename $0` file1[.log] file2[.log] ..."
	echo ""
	echo '	g2psi extracts Standard orientation (or Z-Matrix orientation,'
	echo '	if Standard orientaion does not exist) of molecular coordinates'
	echo '	and MO coefficients from Gaussin output, and then generates'
	echo '	psi88 input files (i.e., psi1, psi2 and psicon).'
	echo '	Molecular coordinates are extracted from the last occurrence'
	echo '	of "Standard orientation" (or "Z-Matrix orientation").'
	echo '	Molecule orbital coefficients are extracted from the last'
	echo '	occurrence of "Molecular Orbital Coefficients".'
}

if [ $# = 0 ]; then printUsage; exit; fi

for FILE in $*
do
	BASE=`basename "$FILE" '\.log'`
	DIR=`dirname "$FILE"`

	if [ ! -r "$FILE" ]
	then
		if [ ! -r "${FILE}.log" ]
		then	echo "File $FILE does not exist."; continue;
		else	FILE=${FILE}.log
		fi
	fi

	nawk '
BEGIN	{
	psi1   = "'$BASE'.psi1";
	psi2   = "'$BASE'.psi2";
	psicon = "'$BASE'.psicon";
	eigen  = 0;
	nbasis = 0;
	basis  = "";
	ndata  = 0;
	natoms = 0;
	nmo    = 0;
	number = "^[-+]?([0-9]+[.]?[0-9]*|[.][0-9]+)([eE][-+]?[0-9]+)?$"
}

/Standard basis:/	{
	basis = $3;
	next;
}

# Standard orientation (or Z-Matrix orientation) is used for overlapping MO.
# Since "Standard orientation" always appears after "Z-Matrix orientation",
# the former is always used if it exists.
/(Standard orientation|Z-Matrix orientation)/	{
	natoms = 0;
	getline; getline; getline; getline;
	while (1) {
		getline;
		if (match($0, "---")) break;
		if ($2 <= 0 + 0) continue;
		natoms++;
		an[natoms] = $2;
		coord[natoms,1] = $3;
		coord[natoms,2] = $4;
		coord[natoms,3] = $5;
	}
	next;
}

# there can be multiple MO coefficients block, e.g., OPT=CALCALL
/Molecular Orbital Coefficients/	{
	nmo++;
	ndata = 0;
	norbs = 0;
	nbasis = 0;
	if (nmo > 1 + 0) {
		close(tmppsi);
		cmd = sprintf("rm %s", tmppsi);
		system(cmd);
	}
	tmppsi = sprintf("%s.tmppsi%d", "'$BASE'", nmo);
	print_psi1_coord(coord, an, natoms, basis, tmppsi);
	next;
}

/EIGENVALUES --/	{
	eigen = 1;
	nf = 0;
# count the number of eigenvalues in this line
	for(i=1;i<=5;i++) {
		eigval = substr($0,(i+1)*10+2,10);
		if (eigval != "          " && eigval != "") nf++;
	}
	norbs += nf;
	next;
}

eigen == 0	{next}

{
# store MO coefficients in a temporary file
	for(i=1;i<=nf;i++) coeff[i] = substr($0,(i+1)*10+2,10) + 0.0;
	if (substr($0, 0, 4) == "    ") {
		for(j=1;j<=nf;j++) {
			for(i=1;i<=nbasis;i++) {
				ndata++;
				printf("% 10.6f", mo[j,i]) > tmppsi
				if (ndata%8==0+0) printf("\n") > tmppsi
			}
		}
		eigen = 0;
		nbasis = 0;
		next;
	}

	nbasis++;
	for(i=1;i<=nf;i++) mo[i,nbasis] = coeff[i];
}

END	{
	if (natoms == 0) {
		printf("\007!!!!! Standard of Z-Matrix orientation is not found. !!!!!\n");
		exit;
	}
#psi1
	printf("\n") > tmppsi
	cmd = sprintf("%s %s %s", "mv", tmppsi, psi1);
	system(cmd);

#psicon
	printf("%s\n", basis) > psicon
	printf(" 1 1 0 1\n")  > psicon
	printf(" 0.075000\n") > psicon

#psi2
	printf("TITLE\n")      > psi2
	printf("SUBTITLE\n")   > psi2
	printf("010000 1.0\n") > psi2
	printf("00\n")         > psi2
	printf("0 COMMENT GOES HERE\n") > psi2
	print_coord(coord, an, natoms, psi2);
	printf("99\n") > psi2
	printf("  61.1000  132.1000    1.1000    0.8500\n") > psi2
	printf("02\n") > psi2
}

function	print_psi1_coord (coord, an, natoms, basis, file,	i)	{
	printf("%s\n", basis)  > file
	printf("AUTO0\n")      > file
	printf("0101  1.0\n")  > file
	printf("0\n")          > file
	print_coord(coord, an, natoms, file)
	printf("99\n") > file
}

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

' $FILE
done