#!/usr/bin/awk -f

# parse OpenMOPAC2009 output for PDB atom labels,
#       produce .mop input marked-up with parenthesized MOPAC comments
# format taken from MOPAC2009, Version 9.007

BEGIN {
 if(!ext) ext="inp"
 
 FS="\n"
# RS="\n *\n *\n *\n *.*\n"
 RS=""

 Nhead=3 # standard MOPAC input header lines (NO extra continued ones considered here!)
 Ndim=3 # dimensions 
}

(FNR==1) {
 if(Natom) printout()
 close(outfile)
 Natom=0
 outfile=FILENAME
 sub(/\.[^.]*$/, "", outfile)
 outfile=outfile "." ext
#print "DEBUG outfile=" outfile
 gotPDB=0
 split("",out) # clear array
 split("",label) # clear array
 
}

(FNR==3) {
	i=2
	while($(++i)!=" *******************************************************************************")
	 continue # skip header intro
#for(j=++i; j<i+Nhead; j++) print "\tDEBUG j=" j "\t'\n" $j  "\n'\n"
	for(j=++i; j<i+Nhead; j++) print $j>outfile
}

/^          CARTESIAN COORDINATES/{
	getline
	getline
	Natom=NF
	for(i=1;i<=Natom;i++) { # split lines into actual fields
		Nsplit=split($i,f," ")
#print "DEBUG 'CARTESIAN COORDINATES' ... Nsplit=" Nsplit "\t$i='\n" $i "\n'\n"
		elem[i]= f[2]
		out[i]= ""
		for(j=1; j<=Ndim; j++) out[i]= out[i] " " sprintf("%7.4f",f[2+j]) " +0 "
#print "DEBUG 'CARTESIAN COORDINATES' ... out[i]=" "\t'\n" out[i] "\n'\n"
	}
}

/^           GEOMETRY AFTER RE-SEQUENCING/{ # special handling...
	getline
	while(NF<3) getline # handle extra blank line 
print "DEBUG 'GEOMETRY AFTER RE-SEQUENCING' ... $0=" "\t'\n" $0 "\n'\n"
		for(i=Nhead+1;i<=NF;i++) { # split lines into actual fields
		if(($i~/^ATOM / || $i~/^HETATM /)){
			Nsplit=split($i,f," ")
			for(j=2; j<=6; j++) label[f[2]]= label[f[2]] " " f[j]
		}else{
			Nsplit=split($i,f,"[)(]")
			elem[i-Nhead]= f[1]
			out[i-Nhead]= f[3]
		}
print "DEBUG 'GEOMETRY AFTER RE-SEQUENCING' ... out[i]=" "\t'\n" out[i] "\n'\n"
	}
}

/^                  Ramachandran Angles/{
print "DEBUG 'Ramachandran Angles' ... #1: NF="NF "\t'\n" $0  "\n'\n"
	getline
	while(NF<3) getline # handle extra blank line preceeding ' COMPUTED CHARGE ON SYSTEM:' 
print "DEBUG 'Ramachandran Angles' ... #2: NF="NF "\t'\n" $0  "\n'\n"
	for(i=1;i<=NF;i++) { # split lines into actual fields
		if(!($i~/^ATOM / || $i~/^HETATM /)) continue
		split($i,f," ")
		for(j=2; j<=6; j++) label[f[2]]= label[f[2]] " " f[j]
print "DEBUG 'Ramachandran Angles' ... label[f[2]]=" "\t'\n" label[f[2]]  "\n'\n"
	}
	gotPDB=1
}

END {
 printout()
}

function printout() {
	for(i=1;i<=Natom;i++){
		out[i]= elem[i] "(" label[i] ") " out[i]
		print out[i] >outfile
	}
	print "\n" >outfile
}

function inchar(i1,i2) {
 Nchar= split($0, char, "")
 string= ""
 for(i=i1; i<=i2; i++) string= string char[i]
 return (string)
}