|
| 1 | +#!/usr/bin/env perl |
| 2 | +################################################################################ |
| 3 | +# |
| 4 | +# Purpose: |
| 5 | +# Reads NWChem's plain stdout (not the ecce_print trace, which carries no |
| 6 | +# thermochemistry and no solvation data) and emits EZEROPT, ETHERM, |
| 7 | +# ENTHALPY, ENTROPY, EGIBBS, THERMOTEMP and, for COSMO, ESOLV. |
| 8 | +# |
| 9 | +# NWChem prints only corrections, so the absolute values ORCA and Gaussian |
| 10 | +# report are rebuilt here. The reference electronic energy is the first |
| 11 | +# "Total <method> energy" after the Hessian banner, i.e. the energy at the |
| 12 | +# analysed geometry, not a later displaced-geometry energy: |
| 13 | +# ETHERM = E + thermal correction to energy |
| 14 | +# ENTHALPY = E + thermal correction to enthalpy |
| 15 | +# ENTROPY = T * S (Hartree, as for ORCA) |
| 16 | +# EGIBBS = ENTHALPY - T*S |
| 17 | +# |
| 18 | +# Temperature = 298.15K |
| 19 | +# Zero-Point correction to Energy = 16.145 kcal/mol ( 0.025729 au) |
| 20 | +# Thermal correction to Energy = 17.922 kcal/mol ( 0.028561 au) |
| 21 | +# Thermal correction to Enthalpy = 18.515 kcal/mol ( 0.029505 au) |
| 22 | +# Total Entropy = 45.039 cal/mol-K |
| 23 | +# |
| 24 | +# ESOLV = sol phase energy - gas phase energy from the COSMO summary. The |
| 25 | +# summary's own "(electrostatic) solvation energy" line is the magnitude |
| 26 | +# of that difference, so the sign is taken from the two totals. |
| 27 | +# |
| 28 | +################################################################################ |
| 29 | + |
| 30 | +$| = 1; |
| 31 | + |
| 32 | +use constant KCAL_PER_HARTREE => 627.5094740631; |
| 33 | + |
| 34 | +my ($e, $temp, $zpe, $dE, $dH, $S, $gas, $sol); |
| 35 | +my $inHessian = 0; |
| 36 | +while (<STDIN>) { |
| 37 | + $inHessian = 1 if (/NWChem Nuclear Hessian and Frequency Analysis/); |
| 38 | + if ($inHessian && !defined($e) |
| 39 | + && /^\s+Total [\w()\-\. ]*? energy\s+=\s+(-?\d+\.\d+)\s*$/) { |
| 40 | + $e = $1; |
| 41 | + } |
| 42 | + if (/^\s*Temperature\s+=\s+(\d+\.\d+)K/) { $temp = $1; } |
| 43 | + elsif (/Zero-Point correction to Energy\s+=.*\(\s*(-?\d+\.\d+)\s+au\)/) { $zpe = $1; } |
| 44 | + elsif (/Thermal correction to Energy\s+=.*\(\s*(-?\d+\.\d+)\s+au\)/) { $dE = $1; } |
| 45 | + elsif (/Thermal correction to Enthalpy\s+=.*\(\s*(-?\d+\.\d+)\s+au\)/) { $dH = $1; } |
| 46 | + elsif (/^\s*Total Entropy\s+=\s+(-?\d+\.\d+)\s+cal\/mol-K/) { $S = $1; } |
| 47 | + elsif (/^\s*gas phase energy\s+=\s+(-?\d+\.\d+)/) { $gas = $1; } |
| 48 | + elsif (/^\s*sol phase energy\s+=\s+(-?\d+\.\d+)/) { $sol = $1; } |
| 49 | +} |
| 50 | + |
| 51 | +sub emit { |
| 52 | + my ($key, $value, $units) = @_; |
| 53 | + print "key: $key\n"; |
| 54 | + print "size:\n1\n"; |
| 55 | + print "values:\n$value\n"; |
| 56 | + print "units:\n$units\n"; |
| 57 | + print "END\n"; |
| 58 | +} |
| 59 | + |
| 60 | +if (defined($e) && defined($dE) && defined($dH) && defined($S) && defined($temp)) { |
| 61 | + my $ts = $S * $temp / 1000.0 / KCAL_PER_HARTREE; |
| 62 | + &emit("EZEROPT", $zpe, "Hartree") if (defined($zpe)); |
| 63 | + &emit("ETHERM", sprintf("%.8f", $e + $dE), "Hartree"); |
| 64 | + &emit("ENTHALPY", sprintf("%.8f", $e + $dH), "Hartree"); |
| 65 | + &emit("ENTROPY", sprintf("%.8f", $ts), "Hartree"); |
| 66 | + &emit("EGIBBS", sprintf("%.8f", $e + $dH - $ts), "Hartree"); |
| 67 | + &emit("THERMOTEMP", $temp, "Kelvin"); |
| 68 | +} |
| 69 | +if (defined($gas) && defined($sol)) { |
| 70 | + &emit("ESOLV", sprintf("%.10f", $sol - $gas), "Hartree"); |
| 71 | +} |
| 72 | +exit(0); |
0 commit comments