#!/usr/bin/perl

# Another fine program brought to you by...
# Michael Sullivan
# sullivan@rsc.anu.edu.au

# "Borrowed" a lot of code from Richard Wong's h298.f

# A program to read in frequencies from a Gaussian log file 
# and calculate ZPVE and Tc with a particular scale factor
# and temperature

# Usage: thermcorr <Gaussian Logfiles>

use strict;
    
foreach my $logfile (@ARGV){
    my @freqs;
    my $pointgroup;
    my $method;
    my $basis;
    my $formula;
    
# Get info from Gaussian Log file    
    ($formula,$pointgroup,$method,$basis,@freqs)=G98LogInfo($logfile);

    my $linear=0;
    my $isLinear=0;
    if ($pointgroup=~/[D|C]\*[H|V]/i){
	$linear=1;
    }
    if ($linear == 1){
	$isLinear ="linear";
    }
    else {
	$isLinear = "not linear";
    }
    
    
#Get user input for Scale factor
    my $scalefactorZPVE=0.8929;
    my $scalefactorTc=$scalefactorZPVE;
    my $scalefactorEntrop=1.0;
    ($scalefactorZPVE,$scalefactorTc,$scalefactorEntrop)=StandScaleFactors($method,$basis);
    my $TEMP = 298.15;
    
    
    print ("This appears to be $formula and is $isLinear.\n");
    print ("The level of theory appears to be $method $basis\n");
    print ("What scale factor should I use for the ZPVE? [$scalefactorZPVE]");
    chomp (my $sfZPVEinput=<STDIN>);
    if ($sfZPVEinput =~ /\d/) {  #* we want to change the scale factor
	$scalefactorZPVE=$sfZPVEinput;
	$scalefactorTc=$scalefactorZPVE;
	$scalefactorEntrop=$scalefactorTc;
    }
    print ("What scale factor should I use for the temperature correction? [$scalefactorTc]");
    chomp (my $sfTcinput=<STDIN>);
#d    print "sfTCinput=$sfTCinput\n";
    if ($sfTcinput =~ /\d/) {  #* we want to change the scale factor
	$scalefactorTc=$sfTcinput;
	$scalefactorEntrop=$scalefactorTc;
    }
#d    print ("scalefactorZPVE=$scalefactorZPVE scalefactorTc=$scalefactorTc\n");
    
    print ("What scale factor should I use for the entropy correction? [$scalefactorEntrop]");
    chomp (my $sfEntropinput=<STDIN>);
    if ($sfEntropinput =~ /\d/) {  #* we want to change the scale factor
	$scalefactorEntrop=$sfEntropinput;
    }
    
# Get the temp
    print ("What temperature do you want this done at? [$TEMP]");
    chomp (my $Temperatureinput=<STDIN>);
    if ($Temperatureinput =~ /\d/) {  #* we want to change the scale factor
	$TEMP=$Temperatureinput;
    }
    
#    undef @scaledfreqs;
    my @scaledfreqs;
    my $numLowModes=0;
    my $numImg=0;
    # Cut off for low modes after scaling
    my $lowmode=260;
    foreach my $indfreq  (@freqs) {
#	print ("$indfreq\n");
#	undef $YN;
	my  $YN;
	my $scaledindfreq=$indfreq*$scalefactorTc;
	if  ($scaledindfreq < 0){
	    print ("Ignoring imaginary mode of $scaledindfreq\n");
	    $YN="Y";
	    $numImg++
	}
	elsif ($scaledindfreq < $lowmode){
	    print ("I found a mode of $scaledindfreq\n");
	    print ("Would you like me to replace this by 1/2 RT? [N] ");
	    chomp ($YN=<STDIN>);
	    if ($YN =~ /^[Yy]/){
		$numLowModes++;
	    }
	}
	unless ($YN =~ /^[Yy]/){
	    push (@scaledfreqs,$scaledindfreq);
	}
    }
#d    $numFreqs=@freqs;
#d    $numSFreqs=@scaledfreqs;
#d    print ("NumFreqs=$numFreqs Num SFreqs=$numSFreqs\n");
    my $PLANCK = 6.6260755E-34;
    my $BOLTZM = 1.380658E-23;
    my $AVOGAD = 6.0221367E23;
    my $SLIGHT = 2.99792458E10;
    my $GASCON = 8.31441;
#    $XKCAL = 4.184;
    
#    print ("$PLANCK    $BOLTZM    $AVOGAD    $SLIGHT    $GASCON    $TEMP      $XKCAL \n");
    # Calculate Zero-point Vibrational Energy
    my $ZPVE=0;
#    print "freqs0=$freqs[0]\n";
    foreach my $indfreq  (@freqs) {
#	print "indfreq=$indfreq\n";
	if  ($indfreq > 0){
	$ZPVE=$ZPVE+($indfreq*$SLIGHT);
	}
    }
    $ZPVE=0.5*$AVOGAD*$PLANCK*$ZPVE/1000;
    my $scaledZPVE=$ZPVE*$scalefactorZPVE;
#d    print ("ZPVE=$ZPVE scaledZPVE=$scaledZPVE\n");
    
# Calculate Temperature Correction
    
    my $SUM = 0.0;
    foreach my $SFREQ  (@scaledfreqs) {
	$SFREQ=$SFREQ*$SLIGHT;
#d	print "SFREQ=$SFREQ ";
	my $XJUNK = exp( $PLANCK * $SFREQ / ($BOLTZM * $TEMP) );
#d	print ("XJUNK=$XJUNK PLANCK=$PLANCK SFREQ=$SFREQ BOLTZM=$BOLTZM TEMP=$TEMP");
	my $TERM = $PLANCK * $SFREQ / ($XJUNK - 1);
#d	print ("TERM =$TERM\n");
	$SUM = $SUM + $TERM;
    }
    my $ROT = $AVOGAD * $SUM;
    my $XFACT;
#d    print ("linear=$linear\n");
    if ($linear == 1) {
	$XFACT = 3.5 ;
    }
    else{
	$XFACT = 4.0  ;
    }
#d    print("XFACT=$XFACT\n");
    my $halfRT=(0.5*$GASCON * $TEMP);
    my $Tc = ($XFACT * $GASCON * $TEMP + $ROT + $numLowModes * $halfRT) / 1000;
#d    print ("RT/2=$halfRT\n");
    
# Calculate Vibrational Entropy
    my $Entropy=0;
#    print "freqs0=$freqs[0]\n";
    foreach my $indfreq  (@freqs) {
#	print "indfreq=$indfreq\n";
	my $SFREQ=$indfreq*$SLIGHT*$scalefactorEntrop;
#d	print "SFREQ=$SFREQ ";
	if ($SFREQ>0){
#	    my $XJUNK = exp( $PLANCK * $SFREQ / ($BOLTZM * $TEMP) );
	    my $PSBT = $PLANCK * $SFREQ / ($BOLTZM * $TEMP);
#	    my $XJUNK = exp;
#	    my $negXJUNK = exp-( $PLANCK * $SFREQ / ($BOLTZM * $TEMP) );
#d	print ("XJUNK=$XJUNK PLANCK=$PLANCK SFREQ=$SFREQ BOLTZM=$BOLTZM TEMP=$TEMP");
#	    my $TERM = $PLANCK * $SFREQ / ($BOLTZM * $TEMP) / (exp($PSBT) - 1) - log(1-exp(-$PSBT));
	    my $TERM = $PSBT / (exp($PSBT) - 1) - log(1-exp(-$PSBT));
#d	print ("TERM =$TERM\n");
	    $Entropy = $Entropy + $TERM *$GASCON;
	}
# S_v= R * Sum[ ( PLANK*SFREQ  /  (BOLTZM*TEMP)  ) /  (XJUNK-1) - ln (1 - 1 / XJUNK) ]
    }

# Print the results
    write;
    format STDOUT =
@<<<<<<<< @####.##K NImag @## Replaced Low Modes @##
$formula, $TEMP,$numImg,$numLowModes
           @#.#### Scaled ZPVE        @####.##### kJ/mol 
$scalefactorZPVE,$scaledZPVE
           @#.#### Scaled Tc          @####.##### kJ/mol
$scalefactorTc,$Tc
           @#.#### Scaled Vib Entropy @####.##### J/mol
$scalefactorEntrop,$Entropy
.
}

sub StandScaleFactors {
    my $method=@_[0];
    my $basis=@_[1];
    my $scalefactorZPVE;
    my $scalefactorTc;
    my $scalefactorEntrop;
    # From Scott and Radom Paper
    if ($method =~ /b3lyp/i) {
	if ($basis =~ /6-31G\(2df,p\)/i){
	    $scalefactorZPVE=0.9854;
	    $scalefactorTc=$scalefactorZPVE;
	    $scalefactorEntrop=$scalefactorZPVE;
	}
	elsif ($basis =~/vtz/i){
	    $scalefactorZPVE=0.985;
	    $scalefactorTc=$scalefactorZPVE;
	    $scalefactorEntrop=$scalefactorZPVE;
	}   
	elsif ($basis =~/6-31G[\(d\)|\*]/i){
	    $scalefactorZPVE=0.9806;
	    $scalefactorTc=0.9989;
	    $scalefactorEntrop=1.0015;
	}   
	else {
	    $scalefactorZPVE=0.985;
	    $scalefactorTc=0.985;
	    $scalefactorEntrop=$scalefactorZPVE;
	}
    }
    elsif ($method =~ /hf/i) {
	if ($basis =~/6-31G[\(d\)|\*]/i){
	    $scalefactorZPVE=0.8929;
	    $scalefactorTc=$scalefactorZPVE;
	    $scalefactorEntrop=$scalefactorZPVE;
	}   
	elsif ($basis =~ /3-21G/i){
	    $scalefactorZPVE=0.9207;
	    $scalefactorTc=0.94444;
	    $scalefactorEntrop=0.9666;
	}
	elsif ($basis =~/6-31\+G[\(d\)|\*]/i){
	    $scalefactorZPVE=0.9153;
	    $scalefactorTc=0.8945;
	    $scalefactorEntrop=0.9027;
	}   
	elsif ($basis =~/6-31G[\(d,p\)|\*\*]/i){
	    $scalefactorZPVE=0.9181;
	    $scalefactorTc=0.8912;
	    $scalefactorEntrop=0.8990;
	}   
	elsif ($basis =~/6-311G[\(d,p\)|\*\*]/i){
	    $scalefactorZPVE=0.9248;
	    $scalefactorTc=0.8951;
	    $scalefactorEntrop=0.9021;
	}   
	elsif ($basis =~/6-311G\(df,p\)/i){
	    $scalefactorZPVE=0.9247;
	    $scalefactorTc=0.8908;
	    $scalefactorEntrop=0.8981;
	}   
	else {
	    $scalefactorZPVE=0.8929;
	    $scalefactorTc=0.8929;
	    $scalefactorEntrop=0.8929;
	}
    }
    # Various 6-31G(d) methods
    elsif ($basis =~ /6-31G[\(d\)|\*]/i){
	if ($method =~/qcisd-fc/i) {
	    $scalefactorZPVE=0.8929;
	    $scalefactorTc=$scalefactorZPVE;
	    $scalefactorEntrop=$scalefactorZPVE;
	}   
	elsif ($method=~/mp2-fu/i){
	    $scalefactorZPVE=0.9661;
	    $scalefactorTc=1.0084;
	    $scalefactorEntrop=1.0228;
	}
	elsif ($method=~/mp2-fc/i){
	    $scalefactorZPVE=0.9670;
	    $scalefactorTc=1.0211;
	    $scalefactorEntrop=1.0444;
	}
	elsif ($method=~/blyp/i){
	    $scalefactorZPVE=1.0126;
	    $scalefactorTc=1.0633;
	    $scalefactorEntrop=1.0670;
	}
	elsif ($method=~/bp86/i){
	    $scalefactorZPVE=1.0108;
	    $scalefactorTc=1.0478;
	    $scalefactorEntrop=1.0527;
	}
	elsif ($method=~/b3p86/i){
	    $scalefactorZPVE=0.9759;
	    $scalefactorTc=0.9864;
	    $scalefactorEntrop=0.9902;
	}
	elsif ($method=~/b3pw91/i){
	    $scalefactorZPVE=0.9774;
	    $scalefactorTc=0.9885;
	    $scalefactorEntrop=0.9920;
	}
	else {
	    $scalefactorZPVE=0.8929;
	    $scalefactorTc=0.8929;
	    $scalefactorEntrop=0.8929;
	}
    }
    else {
	$scalefactorZPVE=0.8929;
	$scalefactorTc=0.8929;
	$scalefactorEntrop=0.8929;
    }
    return($scalefactorZPVE,$scalefactorTc,$scalefactorEntrop);
}
sub G98LogInfo{    
    my $logfile=@_[0];
    my @freqs;
    my $pointgroup;
    my $method;
    my $basis;
    my $formula;
    open(LOG,$logfile) || die ("Cannot open $logfile: $!\n");
    while (<LOG>){
	if (/^ Frequencies --\s+(-?\d*\.\d+)\s*(-?\d*\.\d+)?\s*(-?\d*\.\d+)?\s*/){
#d	    print ("$1\t $2\t$3\n");
	    
	    push (@freqs, $1);
	    if (defined $2) {
		push (@freqs, $2);
	    }
	    if (defined $3) {
		push (@freqs, $3);
	    }
	}
	elsif (/^ Full point group \s+ (\S+)\s+NOp/i){
	    $pointgroup=$1;
#d	    print ("PG=$pointgroup\n");
	}
	elsif (/^ 1\\1\\GINC-\w+\\Freq\\([^\\]*)\\([^\\]*)\\([^\\]*)\\/i){
	    $method=$1;
	    $basis=$2;
	    $formula=$3;
#	    print("method=$method basis=$basis\n");
	}
    }
    if (@freqs eq undef){
	die ("No Frequencies found! Is this a frequency job?\n");
    }
    return($formula,$pointgroup,$method,$basis,@freqs);
}