#!/usr/bin/perl

#Adapted by Dr. Tian Lu (sobereva@sina.com) based on original script
#provided in SI of J. Chem. Theory Comput. 2011, 7, 112–120

use Math::Trig;

### Physical constants ###
$R = 8.314472;		# J K-1 mol-1	Gas constant	
$k = 1.3806504E-23;	# J K-1		Boltzmann constant
$h = 6.62606896E-34;	# J s		Planck constant
$N = 6.02214179E+23;	# mol-1		Avogadro constant
$c = 2.99792458E+10;	# cm s-1	speed of light
$u = 1.660538782E-27;	# kg		amu conversion factor
$a = 101325;		# Pa		atm conversion factor
$e = 2625.49962;	# kJ mol-1	Hartree conversion factor

### Parameters for frequencies ###
$sclz = 0.9770;		# scaling factor for ZPVE
$sclh = 0.9627;		# scaling factor for enthalpy
$scls = 0.9695;		# scaling factor for entropy
$temp = 298.15;		# temperature (K)
$pres = 1;		# pressure (atm)

### HLC and scaling parameters ###
$A1 = 7.1726417223;
$A2 = 7.2642245043;
$B  = 3.6766947822;
$C  = 7.2386727900;
$D  = 2.4044728315;
$E  = 1.0214549544;
$x  = 1.63;
$c1 = 1.3268911659;
$c2 = 0.4032008177;
$c3 = 1.2487668975;
$c4 = 0.4855178788;
$c5 = 1.0772518299;
$c6 = 0.8240421856;

### Help page ###
if (!$ARGV[0]) {
	print ("\nUsage:	g4mp2-6x g09_output_file\(s\)\n\n");
	print ("This script calculates the G4\(MP2\)-6X vibrationless energy, 0 K and 298 K enthalpies \(Hartree\). It checks for normal termination, number of imaginary frequencies, and for near linearity. Gaussian 09 input file should be in the following format:\n");
	print ("
	%chk=h2.chk
	# BMK/6-31+G(2df,p) Opt
	
	Hydrogen molecule G4(MP2)-6X calculation

	0 1
	H        0.000000    0.000000    0.371049 
	H        0.000000    0.000000   -0.371049 

	--Link1--
	%chk=h2.chk
	# Geom=AllCheck Guess=Read BMK/6-31+G(2df,p) Freq

	--Link1--
	%chk=h2.chk
	# Geom=AllCheck Guess=Read CCSD(T,FrzG4)/GTBas1

	--Link1--
	%chk=h2.chk
	# Geom=AllCheck Guess=Read MP2(FrzG4)/GTMP2LargeXP

	--Link1--
	%chk=h2.chk
	# Geom=AllCheck Guess=Read HF/GFHFB3

	--Link1--
	%chk=h2.chk
	# Geom=AllCheck Guess=Read HF/GFHFB4
	\n");

	exit (0);
}

### Script starts here ###
foreach $i (@ARGV) {
	$dead = 0;
	$atyp = 0;
	$name = $i;
	open (LOG, "<$i");
	while (<LOG>) {

	### Determine if the job finished normally ###
		if (/Initial command/) {
			$ter = 0;
		}
		elsif (/Normal termination/) {
			$ter = 1;
		}
		elsif (($ter != 1) && (eof)) {
			$dead = 1;
		}

	### Search for basic molecular properties ###
		elsif (/^ #/) {
			chomp ($Inline = $_);
		}
		elsif (/Charge =.* Multiplicity =/){
			@CMline = split(/\s+/,$_);
		}
		elsif (/Full point group/) {
			@pgline = split(/\s+/,$_);
		}

	### Determine options for G4 frozen core and HLC ###
		elsif ((/^.*Standard orientation:\s*$/) && ($atyp == 0)) {
			$na = 0;
			$nm = 0;
			$nd = 0;
			$nh = 0;
			$ex = 0;
			$line = <LOG>;
			$line = <LOG>;
			$line = <LOG>;
			$line = <LOG>;
			while ($line = <LOG>) {
				last if($line =~ /^\s*------------------+\s*$/);
				@xyz = split(/\s+/,$line);
				$at = $xyz[2];
				if (($at == 11) || ($at == 12) || ($at == 19) || ($at == 20)) {
					$nm++;
				}
				elsif ($at > 20) {
					$nd++;
				}
				elsif ($at == 1) {
					$nh++;
				}
				elsif (($at == 2) ||($at == 5) || ($at == 13)) {
					$ex = 1;
				}
				$na++;
			}
			$atyp = 1;
		}

	### Search for frequencies and apply scale factors ###
		elsif (/Harmonic frequencies/) {
			$vib = 0;
			$nvib = 0;
			$imag = 0;
			$pg = $pgline[4];
			$mult = @CMline[6];
			if ($pg =~ /OH/) {
				$moltype = 0;	# atom
			}
			elsif ($pg =~ /[D|C]\*[H|V]/) {
				$moltype = 1;	# linear molecule
			}
			else {
				$moltype = 2;	# non-linear molecule
			}
		}
		elsif (/Frequencies/) {
			@freqs = split(/\s+/,$_);
			if (@freqs[3] > 0) {
				$freq{$vib} = @freqs[3];
				$zfreq{$vib} = $freq{$vib} * $sclz;
				$hfreq{$vib} = $freq{$vib} * $sclh;
				$sfreq{$vib} = $freq{$vib} * $scls;
				$vib++;
				$nvib++;
			}
			elsif (@freqs[3] < 0) {
				$imag++;
			}
			if (@freqs[4] > 0) {
				$freq{$vib} = @freqs[4];
				$zfreq{$vib} = $freq{$vib} * $sclz;
				$hfreq{$vib} = $freq{$vib} * $sclh;
				$sfreq{$vib} = $freq{$vib} * $scls;
				$vib++;
				$nvib++;
			}
			elsif (@freqs[4] < 0) {
				$imag++;
			}
			if (@freqs[5] > 0) {
				$freq{$vib} = @freqs[5];
				$zfreq{$vib} = $freq{$vib} * $sclz;
				$hfreq{$vib} = $freq{$vib} * $sclh;
				$sfreq{$vib} = $freq{$vib} * $scls;
				$vib++;
				$nvib++;
			}
			elsif (@freqs[5] < 0) {
				$imag++;
			}
		}

	### Search for other physical properties ###
		elsif (/Molecular mass/) {
			@massline = split(/\s+/,$_);
			$mass = $massline[3];
		}
		elsif (/Rotational symmetry number/) {
			@rotsymline = split (/\s+/,$_);
			$rotsym = $rotsymline[4];
		}
		elsif (/Rotational temp/) {
			@rottempline = split (/\s+/,$_);
			$rottempx = $rottempline[4];
			$rottempy = $rottempline[5];
			$rottempz = $rottempline[6];
			if ($rottempx =~ /\*/ || $rottempy =~ /\*/ || $rottempz =~/\*/) {
				$moltype = 3;	# supposely linear molecule
			}
		}

	### Search for CCSD(T)/GTBas1 component energies and number of electrons ###
		
		elsif ((/SCF Done/) && ($Inline =~ /CCSD.*GTBas1/i)) {
			@hf = split(/\s+/,$_);
			$hf_bs1 = $hf[5];
		}
		elsif ((/NOA=.*NOB=/) && ($Inline =~ /CCSD.*GTBas1/i)) {
			@noe = split (/\s+/,$_);
			$noa = $noe[4] - $nm * 4 - $nd * 5;
			$nob = $noe[6] - $nm * 4 - $nd * 5;
		}
		elsif ((/alpha-alpha/) && ($Inline =~ /CCSD.*GTBas1/i)) {
			@aa = split (/\s+/,$_);
			$aa[6] =~ s/D/E/g;
			$aa_bs1 = $aa[6];
		}
		elsif ((/alpha-beta/) && ($Inline =~ /CCSD.*GTBas1/i)) {
			@ab = split (/\s+/,$_);
			$ab[6] =~ s/D/E/g;
			$ab_bs1 = $ab[6];
		}
		elsif ((/beta-beta/) && ($Inline =~ /CCSD.*GTBas1/i)) {
			@bb = split (/\s+/,$_);
			$bb[6] =~ s/D/E/g;
			$bb_bs1 = $bb[6];
		}
		elsif ((/E\(CORR\)=/) && ($Inline =~ /CCSD.*GTBas1/i)) {
			@cc = split (/\s+/,$_);
			$cc[4] =~ s/D/E/g;
			$cc_bs1 = $cc[4] - $hf_bs1;
		}
		elsif ((/CCSD\(T\)= /) && ($Inline =~ /CCSD.*GTBas1/i)) {
			@t = split (/\s+/,$_);
			$t[2] =~ s/D/E/g;
			$te_bs1 = $t[2] - ($cc_bs1 + $hf_bs1);
		}

	### Search for MP2/GTMP2LargeXP component energies ###
		elsif ((/SCF Done/) && ($Inline =~ /MP2.*GTMP2LargeXP/i)) {
			@hf = split(/\s+/,$_);
			$hf_lxp = $hf[5];
		}
		elsif ((/alpha-alpha/) && ($Inline =~ /MP2.*GTMP2LargeXP/i)) {
			@aa = split (/\s+/,$_);
			$aa[6] =~ s/D/E/g;
			$aa_lxp = $aa[6];
		}
		elsif ((/alpha-beta/) && ($Inline =~ /MP2.*GTMP2LargeXP/i)) {
			@ab = split (/\s+/,$_);
			$ab[6] =~ s/D/E/g;
			$ab_lxp = $ab[6];
		}
		elsif ((/beta-beta/) && ($Inline =~ /MP2.*GTMP2LargeXP/i)) {
			@bb = split (/\s+/,$_);
			$bb[6] =~ s/D/E/g;
			$bb_lxp = $bb[6];
		}

	### Search for HF/GFHFB3 energy ###
		elsif ((/SCF Done/) && ($Inline =~ /HF.*GFHFB3/i)) {
			@hf = split(/\s+/,$_);
			$hf_hfb3 = $hf[5];
		}

	### Search for HF/GFHFB4 energy ###
		elsif ((/SCF Done/) && ($Inline =~ /HF.*GFHFB4/i)) {
			@hf = split(/\s+/,$_);
			$hf_hfb4 = $hf[5];
		}       
	}

	printf ("%-40s ", "$name");

	if (($dead != 1) && (eof)) {

	### Calculate themochemistry ###
		### ZPVE ###
		$sumzfreq = 0;
		for ($j = 0; $j < $nvib; $j++) {
			$sumzfreq = $sumzfreq + $zfreq{$j};
		}
		$zpve = (1/2 * $h * $c * $N * $sumzfreq / 1000) / $e;

		### Translational ###
		$qt = (2 * pi * $mass * $u * $k * $temp / $h**2)**(3/2) * $k * $temp / ($pres * $a);
		$St = $R * (log($qt) + 1 + 3/2);
		$Et = 3/2 * $R * $temp;
	
		### Electronic ###
		$Se = $R * log($mult);
	
		### Rotational ###
		if ($moltype == 0) {
			$Sr = 0;
			$Er = 0;
		}
		elsif ($moltype == 1) {
			$qr = 1 / $rotsym * $temp / $rottempx;
			$Sr = $R * (log($qr) + 1);
			$Er = $R * $temp;
		}
		elsif ($moltype == 2) {
			$qr = pi**(1/2) / $rotsym * ($temp**(3/2) / ($rottempx * $rottempy * $rottempz)**(1/2));
			$Sr = $R * (log($qr) + 3/2);
			$Er = 3/2 * $R * $temp;
		}

		### Vibrational ###
		$Sv = 0;
		for ($j = 0; $j < $nvib; $j++) {
			$Qvs{$j} = $h * $sfreq{$j} * $c / $k;
			$Sv = $Sv + ($Qvs{$j} / $temp) / (exp($Qvs{$j} / $temp) - 1) - log(1 - exp(-1 * $Qvs{$j} / $temp));
		}
		$Sv = $R * $Sv;

		$Ev = 0;
		for ($j = 0; $j < $nvib; $j++) {
			$Qvh{$j} = $h * $hfreq{$j} * $c / $k;
			$Ev = $Ev + $Qvh{$j}  / (exp($Qvh{$j} / $temp) - 1);
		}
		$Ev = $R * $Ev;

		$H = (($Et + $Er + $Ev + $R * $temp) / 1000) / $e;
		$Stot = ($St + $Se + $Sr + $Sv) / $e; 

	### Calculate G4(MP2)-6X electronic energy ###
		### Total additivity energy ###
		$hf_cbs = ($hf_hfb4 - $hf_hfb3 * exp(-$x)) / ( 1 - exp(-$x));
		$scse2_1 = $c1 * $ab_bs1 + $c2 * ($aa_bs1 + $bb_bs1);
		$scse2_2 = $c3 * $ab_lxp + $c4 * ($aa_lxp + $bb_lxp);
		$scc = $c5 * $cc_bs1 - $scse2_1;
		$st = $c6 * $te_bs1;
		$g4a = $hf_cbs + $scse2_2 + $scc + $st;

		### HLC for molecules ###
		if ($na > 1) {
			if ($mult == 1) {
				if (($nob != 1) || (($nob == 1) && ($nh != 0))) {
					$hlc = -$A1 * $nob;
				}
				elsif (($nob == 1) && ($nh == 0)) {
					$hlc = -$E *$nob;
				}
			}
			else {
				$hlc = -$A2 * $nob - $B * ($noa - $nob);
			}
		}

		### HLC for atomic species ###
		elsif ($na == 1) {
			if ($mult == 1) {
				if (($nob != 1) || (($nob == 1) && ($ex != 0))) {
					$hlc = -$C * $nob;
				}
				elsif (($nob == 1) && ($ex == 0)) {
					$hlc = -$E *$nob;
				}
			}
			else {
				 $hlc = -$C * $nob - $D * ($noa - $nob);
			}
		}

		$hlc = $hlc / 1000;
		$g4e = $g4a + $hlc;
		$g40 = $g4e + $zpve;
		$g4298 = $g40 + $H;
		$U = $g4298 - $temp*$R/1000/$e;
		$Gibbs = $g4298 - $Stot * $temp/1000;

	### Print results ###
		printf ("\n%18s %8.2f", "Temperature (K)", "$temp");
		if ($moltype != 3) {
			printf ("\n%8s %5d", "NImag", "$imag");
			printf ("\n%8s %16.8f", "E_ele", "$g4e");
			printf ("\n%8s %16.8f", "H(0K)", "$g40");
			printf ("\n%8s %16.8f", "H(298K)", "$g4298");
			printf ("\n%8s %16.8f", "U(298K)", "$U");
			printf ("\n%8s %16.8f", "G(298K)", "$Gibbs");
		}
		else {
			print ("		    ------ Is this linear? ------");
		}
	}

	elsif (($dead == 1) && (eof)) {
			print ("			  Error termination");	
	}

	print ("\n");
	close (LOG);
}
print ("\n");