#!/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 (
) { ### 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 =
; $line =
; $line =
; $line =
; while ($line =
) { 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");