#!/usr/local/bin/perl

###############################################################

# Program to read a jaguar restart file and output a molden molf file
#
# Version History
# 12/3/98  Version 1.00b    Beta release
#
#
# Bug List
# (1) Need to find the jaguar basis set file more intelligently
# (2) Need to use the backup record to find the basis set rather than just
#     copping out and using 6-31.
# (3) Need to support z-matrices

use Getopt::Std;

# This line may have to be changed to reflect the local configuration:
$basis_file_name = '/source/schrodinger/v35msc/jaguar/data/default.basis';

if (!@ARGV) {
    print "Usage: jag2molf [options] file.in\n\n";
    print "Options:\n";
    print "-c X   Set the wave function cutoff to X (default 0.01)\n";
    print "-h     Print the help message\n";
    print "-v     Verbose Mode\n";
    die "\n";
}

getopts('c:hv');

if ($opt_h) {
    print "Usage: jag2molf [options] file.in\n\n";
    print "Options:\n";
    print "-c X   Set the wave function cutoff to X (default 0.01)\n";
    print "-h     Print the help message\n";
    print "-v     Verbose Mode\n";
    die "\n";
}

$input_file = shift(@ARGV);

if ($opt_v) {
    $verbose = 1;
} else {
    $verbose = 0;
}

if ($opt_c) {
    $cutoff = $opt_c;
} else {
    $cutoff = 0.01;
}


$output_file = $input_file;
$output_file =~ s/\..*//;
$output_file = $output_file . ".molf";

open(IN,"< $input_file") || print "Can\'t open file $input_file\n";

get_atomic_parameters();

$got_orbs = 0;
$got_geo = 0;
$got_basis = 0;

while(<IN>){

    if (/[\$\&]zmat/) {
	get_geo();
    } elsif (/guess/) {
	($basis_name) = /basgss=(\S+)/;

	$ipolar = set_polar($basis_name);
	$idiffuse = set_diffuse($basis_name);

	$basis_name = clean_basis_name($basis_name);

	get_guess();
    }
}

get_basis();

write_molden();

###############################################################

sub set_polar{

    $bname = shift;

    $_ = $bname;
    if (/\*\*/) {
	if ($verbose) {
	    print "Polarization functions on all atoms\n";
	}
	$ipol = 2;
    } elsif (/\*/) {
	if ($verbose) {
	    print "Polarization functions only on heavy atoms\n";
	}
	$ipol = 1;
    } else {
	if ($verbose) {
	    print "No polarization functions\n";
	}
	$ipol = 0;
    }
    
    return $ipol;    

}

###############################################################

sub set_diffuse{

    $bname = shift;

    $_ = $bname;
    if (/\+\+/) {
	if ($verbose) {
	    print "Diffuse functions on all atoms\n";
	}
	$idiff = 2;
    } elsif (/\+/) {
	if ($verbose) {
	    print "Diffuse functions only on heavy atoms\n";
	}
	$idiff = 1;
    } else {
	if ($verbose) {
	    print "No diffuse functions\n";
	}
	$idiff = 0;
    }
    
    return idiff;

}


###############################################################

sub clean_basis_name{

    $bname = shift;

    $bname = uc $bname;

    $bname =~ s/G//; #remove the G's from e.g. 3-21G
    $bname =~ s/\++//; #remove the +, ++
    $bname =~ s/\*+//; #remove the *, **
    
    return $bname;
}

###############################################################

sub get_geo{

# Only works for Cartesian now

    $nat = 0;
    while (<IN>) {
	last if (/\&/);
	last if (/\$/);
	($atom[$nat],$x[$nat],$y[$nat],$z[$nat]) = split;
	$nat++;

    }    

    $got_geo = 1;

}

###############################################################

sub get_guess{

    $norb = -1;
    LINE: while(<IN>) {
	last LINE if (/\$/);
	last LINE if (/\&/);
	if (/Orbital Energy/) { # Get the orbital header info
	    $norb++;
	    $nbf = 0;
	    ($orb_energy[$norb]) = /\s+(\-*\d*\.\d*)\s+Occupation/;
	    ($orb_occ[$norb]) = /Occupation\s+(\d*\.\d*)/;
	    ($sym_label[$norb]) = /Symmetry\s+(\w+)/;
	} else {
	    ($orb[$norb][$nbf],$orb[$norb][$nbf+1],$orb[$norb][$nbf+2],
	     $orb[$norb][$nbf+3]) = split;
	    
	    if ($orb[$norb][$nbf+1] eq ''){
		$nbf++;
	    } elsif ($orb[$norb][$nbf+2] eq '') {
		$nbf += 2;
	    } elsif ($orb[$norb][$nbf+3] eq '') {
		$nbf += 3;
	    } else {
		$nbf += 4;
	    }
	}

    }

    $norb++; # caused by the stupid way we're incrementing
    $got_orbs = 1;
	    
}

###############################################################

sub get_basis{

    $nbf_mo = $nbf;

    $nbf = 0;

    for ($iat = 0; $iat < $nat; $iat++) {

	$element = uc $atom[$iat]; #everything uppercase
	$element =~ s/\W.*//;   
	$element =~ s/\d.*//; 

	$istart[$iat] = $nbf;

	$worked = get_atom_basis();

# Desperate attempt to get a backup function
	if (!$worked) {
	    get_631_basis();
	}

	$iend[$iat] = $nbf;

    }

    $got_basis = 1;

}

###############################################################

sub get_atom_basis{

    open(BASIS,"<$basis_file_name") || die "Can\'t open $basis_file_name\n";
    $igotbf = 0;
    FETCH: while(<BASIS>) {
	if (/BASIS/) {
	    @basis_line = split/,|\s+/;

	    $basis_match = parse_bline($basis_name);

	    if ($basis_match) {


		if ($verbose) {
		    print "BASIS MATCH FOUND!\n";
		}

		while(<BASIS>) {
		
		    @el_line = split;
		    $el_match = parse_eline();

		    if ($el_match) {
			read_basis_element();
			close(BASIS);
			last FETCH;
		    } else {
			skip_basis_element();
		    }
		}
	    }  # if ($basis_match)

	} # if (/BASIS/)
    } # while(<BASIS>)

    return $igotbf;
}

##############################################################

sub get_631_basis{

    open(BASIS,"<$basis_file_name") || die "Can\'t open $basis_file_name\n";
    $igotbf = 0;
    FETCH: while(<BASIS>) {
	if (/BASIS/) {
	    @basis_line = split/,|\s+/;

	    $basis_match = parse_bline('6-31');

	    if ($basis_match) {


		if ($verbose) {
		    print "BASIS MATCH FOUND!\n";
		}

		while(<BASIS>) {
		
		    @el_line = split;
		    $el_match = parse_eline();

		    if ($el_match) {
			read_basis_element();
			close(BASIS);
			last FETCH;
		    } else {
			skip_basis_element();
		    }
		}
	    }  # if ($basis_match)

	} # if (/BASIS/)
    } # while(<BASIS>)

    return $igotbf;
}

##############################################################

sub parse_eline{

    $nentries = scalar(@el_line);
    if ($nentries == 1) {
	$entry = uc $el_line[0];
	if ($verbose){
	    print " Checking if $entry matches $element...";
	}
	if ($entry eq $element) {
	    if ($verbose) {
		print "YES!\n";
	    }
	    return 1;
	} else {
	    if ($verbose) {
		print "NO.\n";
	    }
	    return 0;
	}
    } else {
	return 0;
    }

    return 0; # just in case...
}

##############################################################

sub parse_bline{

    my $bn = shift;

# This routine parases the BASIS line to figure out 
#   (1) whether the line matches the basis set
#   (2) ECP, backup, etc.

    $got_term = 0;
    $ecp = 0;

    $got_backup = 0;
    $num_backup = 0;
    $num_bnames = 0;

    foreach $entry (@basis_line) {

	$entry = clean_basis_name($entry);

	if ($entry eq 'BASIS') {
	    # Do nothing here
	} elsif ($entry eq '5D' or $entry eq '6D') {
	    $got_term = 1;
	} elsif (!$got_term) {
	    $bnames[$num_bnames] = $entry;
	    $num_bnames++;
	} elsif ($entry eq 'ECP') {
	    $ecp = 1;
	} elsif ($entry eq 'BACKUP'){
	    $got_backup = 1;
	} elsif ($got_backup) {
	    $backup[$num_backup] = $entry;
	    $num_backup++;
	} else {
	    print "WARNING: don\'t know what to do with $entry ";
	    print "in parse_bline\n";
	}

    }

    foreach $bname (@bnames) {
	if ($verbose) {
	    print "Does $bname match $bn?     ";
	}
	if ($bname eq $bn) {
	    if ($verbose) {print "YES!\n"};
	    return 1;
	}
	if ($verbose) {print "NO.\n"};
    }

    return 0;

}

##################################################################

sub read_basis_element{

    if ($verbose) {
	print "Reading element $element in basis set $basis_name\n";
    }

    while(<BASIS>) {
	last if (/\*\*/);
	($shell_type[$nbf], $func_type,$bs1,$bs2,$bs3) = split;
	
	if ($bs3 eq "-") {
	    $nprim[$nbf] = $bs1 + $bs2;
	} elsif ($bs2 eq "-") {
	    $nprim[$nbf] = $bs1;
	} else {
	    die "Error reading basis file, line = \n$_\n";
	}

	if ($verbose) {
	    print "$nbf $shell_type[$nbf] $func_type $nprim[$nbf]\n";
	}

	$igotbf=1;
			
			
	if ($func_type == 0) {  #normal basis function
	    read_prims();
	} elsif ($func_type > 0) { #polarization function
	    if ($ipolar >= $func_type) {
		read_prims();
	    } else {
		skip_prims();
	    }
	} elsif ($func_type < 0) { #diffuse function
	    $ft_pos = abs $func_type;
	    if ($idiffuse  >= $ft_pos) {
		read_prims();
	    } else {
		skip_prims();
	    }
	} else {
	    die "Bad function type $func_type\n";
	}

    };

}

##############################################################

sub skip_basis_element{

    $igotbf = 0;
    while (<BASIS>) {
	last if (/\*\*/);
	last if (/BASIS/);
    };
    
}

###############################################################

sub read_prims{

    for ($i = 0; $i<$nprim[$nbf]; $i++) {
	$_ = <BASIS>;

	if ($shell_type[$nbf] eq 'SP') {
	    ($exp[$nbf][$i], 
	     $cont_coef[$nbf][$i], 
	     $cont_coefp[$nbf][$i]) = split;
	} else {
	    ($exp[$nbf][$i], 
	     $cont_coef[$nbf][$i]) = split;
	}
	if ($verbose) {
	    print "    $exp[$nbf][$i] $cont_coef[$nbf][$i] ";
	    print " $cont_coefp[$nbf][$i]\n";
	}
    }
    $nbf++;
    $igotbf = 1;

}

###############################################################

sub skip_prims{
    for ($i = 0; $i<$nprim[$nbf]; $i++) {
	$_ = <BASIS>;
    }
}


###############################################################

sub write_molden{

    open(OUT,">$output_file") || die "Can\'t open $output_file\n";


    if ($got_geo) { 
	write_molden_atoms();
    }

    if ($got_basis) {
	write_molden_basis();
    }

    if ($got_orbs) {
	write_molden_orbs();
    }
}

###############################################################

sub write_molden_atoms{

# Atoms Section
    print OUT "[Atoms] Angs\n";
    for ($iat = 0; $iat < $nat; $iat++) {

	$element = uc $atom[$iat]; #everything uppercase
	$element =~ s/\W.*//;   
	$element =~ s/\d.*//; 

	$iatp1 = $iat+1;

	printf OUT "%s %5.0f %5.0f %10.4f %10.4f %10.4f\n",
	$element, $iatp1, $atno{$element},$x[$iat],$y[$iat],$z[$iat];

    }

}

###############################################################

sub write_molden_basis{

# Basis Set Section
    print OUT "[GTO]\n";
    for ($iat = 0; $iat < $nat; $iat++) {

	$iatp1 = $iat + 1;

	print OUT "$iatp1 0\n";
	for ($ibf = $istart[$iat]; $ibf < $iend[$iat]; $ibf++) {

	    if ($shell_type[$ibf] eq 'SP') {
		print OUT "S   $nprim[$ibf]  1.00\n";
		for ($iprim = 0; $iprim < $nprim[$ibf]; $iprim++){
		    printf OUT " %10.4f %10.4f  \n", 
		    $exp[$ibf][$iprim], $cont_coef[$ibf][$iprim];
		}
		print OUT "P   $nprim[$ibf]  1.00\n";
		for ($iprim = 0; $iprim < $nprim[$ibf]; $iprim++){
		    printf OUT " %10.4f %10.4f  \n", 
		    $exp[$ibf][$iprim], $cont_coefp[$ibf][$iprim];
		}
	    } else {
		print OUT "$shell_type[$ibf]  $nprim[$ibf] 1.00\n";
		for ($iprim = 0; $iprim < $nprim[$ibf]; $iprim++) {
		    printf OUT " %10.4f %10.4f  \n", 
		    $exp[$ibf][$iprim], $cont_coef[$ibf][$iprim];
		}
	    }
	}
	print OUT "\n";
    }
    print OUT "\n";

}

###############################################################

sub write_molden_orbs{

# Orbital Section
    print OUT "[MO] \n";
    for ($iorb = 0; $iorb < $norb; $iorb++) {

	print OUT "Sym=$sym_label[$iorb]\n";
	print OUT "Ene=$orb_energy[$iorb]\n";
	print OUT "Spin=Alpha\n";
	print OUT "Occup=$orb_occ[$iorb]\n";

	for ($ibf = 0; $ibf < $nbf_mo; $ibf++) {
	    if (abs($orb[$iorb][$ibf]) >= $cutoff) {
		$ibfp1 = $ibf + 1;
		printf OUT "%5.0f %10.4f\n",$ibfp1,$orb[$iorb][$ibf];
	    }
	}
    }

}


###############################################################

sub get_atomic_parameters{

    %atno = (
	     H => 1,
	     HE => 2,
	     LI => 3,
	     BE => 4,
	     B => 5,
	     C => 6,
	     N => 7,
	     O => 8,
	     F => 9,
	     NE => 10,
	     NA => 11,
	     MG => 12,
	     AL => 13,
	     SI => 14,
	     P => 15,
	     S => 16,
	     CL => 17,
	     AR => 18,
	     K => 19,
	     CA => 20,
	     SC => 21,
	     TI => 22,
	     V => 23,
	     CR => 24,
	     MN => 25,
	     FE => 26,
	     CO => 27,
	     NI => 28,
	     CU => 29,
	     ZN => 30,
	     GA => 31,
	     GE => 32,
	     AS => 33,
	     SE => 34,
	     BR => 35,
	     KR => 36,
	     RB => 37,
	     SR => 38,
	     Y => 39,
	     ZR => 40,
	     NB => 41,
	     MO => 42,
	     TC => 43,
	     RU => 44,
	     RH => 45,
	     PD => 46,
	     AG => 47,
	     CD => 48,
	     IN => 49,
	     SN => 50,
	     SB => 51,
	     TE => 52,
	     I => 53,
	     XE => 54,
	     CS => 55,
	     BA => 56,
	     LA => 57,
	     CE => 58,
	     PR => 59,
	     ND => 60,
	     PM => 61,
	     SM => 62,
	     EU => 63,
	     GD => 64,
	     TB => 65,
	     DY => 66,
	     HO => 67,
	     ER => 68,
	     TM => 69,
	     YB => 70,
	     LU => 71,
	     HF => 72,
	     TA => 73,
	     W => 74,
	     RE => 75,
	     OS => 76,
	     IR => 77,
	     PT => 78,
	     AU => 79,
	     HG => 80,
	     TL => 81,
	     PB => 82,
	     BI => 83,
	     PO => 84,
	     AT => 85,
	     RN => 86,
	     FR => 87,
	     RA => 88,
	     AC => 89,
	     TH => 90,
	     PA => 91,
	     U => 92,
	     NP => 93,
	     PU => 94,
	     AM => 95,
	     CM => 96,
	     BK => 97,
	     CF => 98,
	     ES => 99,
	     FM => 100,
	     MD => 101,
	     NO => 102,
	     LR => 103 );
}

###############################################################