Schupf Computational Chemistry Lab

#!/usr/bin/perl # format DeFt and GAMESS output file in HTML $filert = shift(@ARGV) ; # Extract first argument into pattern $file = $filert.".dft" ; # DeFT output file open(DFT,$file) || die "Can't find file $file\n"; $file = $filert.".gms" ; # GAMESS output file open(GMS,$file) || die "Can't find file $file\n"; # get angles from dft z-matrix while ( ) { last if /esp/ ; if ( /total energy/ ) { chop($E = $_) ; $E =~ s/^.*=// ; $E =~ s/ *// ; } # end if total energy for cycle next unless /z-matrix/ ; $j = 0 ; $_ = ; $_ = ; $_ = ; while ( ) { last if /^\n/ ; $j++ ; if ( $j == 1 ) { ($x,$el,$r,$angle[$j],$ael[$j],$ael3[$j],$x) = split } else { ($x,$x,$x,$angle[$j],$x,$ael[$j],$ael3[$j],$x) = split } } # end while z-matrix } # end while optimization steps $nang = $j ; # get charge $_ = ; $j = 0 ; while ( ) { last unless /esp/ ; $j++ ; ($x,$q[$j]) = split("=") ; $q[$j] = substr($q[$j],1,6) ; } # end while charge # get mo orbital energies while ( ) { last if /kohn/ } $_ = ; # skip blank line # get alpha orbital energies $j = 0 ; while ( ) { last if /^\n/ ; $j++ ; ($x,$mo[$j],$x,$moe[$j],$x,$occa[$j]) = split('[#,=]') ; chop($moe[$j]) ; chop($moe[$j]) ; # drop eV units $moe[$j] =~ s/ *//g ; $homo = $j if $occa[$j] == 1 ; } # end while alpha orbitals $imo = $j ; # get beta orbital occupations $j = 0 ; while ( ) { last if /^\n/ ; $j++ ; ($x,$x,$x,$x,$x,$occb[$j]) = split('[#,=]') ; } # end while alpha orbitals for($i=$j+1;$i<=$imo;$i++) { $occb[$i] = 0 } # get dipole moment while ( ) { last if /total :/ } ($x,$x,$x,$x,$dpm) = split ; chop($dpm) ; # find multiplicity and charge if not default $mult = 1 ; $charge = 0 ; while ( ) { last if /ECHO/ } while ( ) { last unless /INPUT/ ; if (/MULT/) { /MULT=[ ]*([^ ]*)/; $mult = $1 ; } if (/ICHARG/) { /ICHARG=[ ]*([^ ]*)\b/; $charge = $1 ; } } # end unless specifications # get element names while ( ) { last if /CHARGE/ } $iel = 0 ; while ( ) { last if /^\n/ ; $iel++ ; ($x,$el[$iel],$x) = split ; } # end while elements # get bond lengths and bond orders bond: while ( ) { last if /BOND ORDER/ } while ( ) { last if /ATOM/ } $ib = 0 ; while ( ) { last if /^\n/ ; $ib++ ; ($x,$a1[$ib],$a2[$ib],$d[$ib],$bo[$ib],$a1[$ib+1],$a2[$ib+1],$d[$ib+1], $bo[$ib+1],$a1[$ib+2],$a2[$ib+2],$d[$ib+2],$bo[$ib+2]) = split ; if ( $a1[$ib+1] eq '' ) { last } if ( $a1[$ib+2] eq '' ) { $ib++ ; last } $ib += 2 ; } # end unless bonds # #print "Content-type: text/html\n\n" ; print "\n\nSchupf Computational Chemistry Lab\n"; print "\n\n"; print "\n" ; # get structures for display $file = $filert.".htm" ; # html output file with stick and embedded mol src open(HTM,$file) || die "Can't find file $file\n"; while ( ) { print } # print "The ion charge is $charge. " if $charge ; print"The multiplicity is $mult.
\n" if $mult != 1 ; print "

\n" if $charge != 0 || $mult != 1; print "Tell me about the atomic charges, dipole moment,\n" ; print "bond lengths, \n" ; print "angles, \n" ; print "bond orders,
\n" ; print "molecular orbital energies,\n" ; print "or total energy.
\n" ; #print " or hybridization.\n" ; print "Tell me about the best Lewis structure.\n"; print "\n" ; print "

Atomic Charges and Dipole Moment

\n"; for ($k=1;$k<=$iel;$k++) { $nm1 = $el[$k] ; print "$nm1$k charge=$q[$k]
\n"; } # continue charge $dpm = 0 if $dpm < 0.0001 ; print " with a dipole moment of $dpm Debye\n" ; print "

\n"; print "\n" ; print "

Bond Lengths:

\n"; for ($k=1;$k<=$ib;$k++) { $nm1 = $el[$a1[$k]] ; $nm2 = $el[$a2[$k]] ; print "between $nm1$a1[$k] and $nm2$a2[$k]: distance=$d[$k] ang___\n"; print "
\n" unless $k%2 ; } # continue bonds print "

\n"; if ( $nang > 0 ) { print "\n" ; print "

Bond Angles:

\n"; for ($k=1;$k<=$nang;$k++) { $ae1 = $k+2 ; $x = substr($angle[$k],0,5) ; $nm1 = $el[$ae1] ; $nm2 = $el[$ael[$k]] ; $nm3 = $el[$ael3[$k]] ; print "for $nm1$ae1-$nm2$ael[$k]-$nm3$ael3[$k]: angle=$x deg___\n"; print "
\n" unless $k%2 ; } # continue bonds print "

\n"; print "Top of page.\n" ; } # end if angles print "\n" ; print "

Bond Orders (Mulliken):

\n"; for ($k=1;$k<=$ib;$k++) { $nm1 = $el[$a1[$k]] ; $nm2 = $el[$a2[$k]] ; print "between $nm1$a1[$k] and $nm2$a2[$k]: order=$bo[$k]___\n"; print "
\n" unless $k%2 ; } # continue bonds print "

\n"; print "Top of page.\n" ; # get hybridization for each bond print "\n" ; print "

Best Lewis Structure

\n"; print "The Lewis structure that is closest to your structure is determined.\n"; print "The hybridization of the atoms in this idealized Lewis structure \n"; print "is given in the table below. \n"; # check for resonance structures while ( ) { last if /Thresh/ } $_ = ; while ( ) { last if /----/ ; ($x,$j,$k,$l) = split ; } # end while resonance $deloc = 0 ; # true false for delocaized structure if ( ( $k < 1.9 && $mult ==1 ) || ( $k < 0.90 ) ) { $deloc = 1 } while ( ) { last if /----/ ; $deloc = 1 if /delocalized structure/ ; } # end while delocalized if ( $deloc ) { print "Please note that your structure can't be well described by a single\n"; print " Lewis structure, because of extensive delocalization.\n

\n"; } # end if delocalized $updn = 1 ; $updn = 2 if $mult > 1 ; if ( $updn > 1 ) { print "The Lewis structure is built for the up and down electrons, \n"; print "separately. Note that the up and down structures can be very \n"; print "different.
\n" ; } # end if up and down explain print "

\n" ; print "$j A $type orbital for $nm1 with $occ electrons
\n"; $bond = 1 ; # for print statement below $qual = -2 ; # need two more lines of input next ; } # end if bond header # for lone pair orbitals if ( /LP/ ) { ($j,$occ,$type,$k,$el,$s,$i,$p,$i,$d) = split(/[()]/) ; next if $occ < 0.10 ; chop($occ) ; $el =~ s/ *//g ; chop($el) ; $j =~ s/ *// ; chop($s) ; chop($p) ; chop($d) ; print "

\n" ; print "$j A lone pair orbital for $el with $occ electrons
\n"; $bond = 0 ; #for print statement below $qual = 0 ; # don't need any more lines of input } # end if LP header else { next unless $qual ; $qual++ ; ($x,$pol,$k,$el,$s,$j,$p,$i,$d) = split(/[()\*]/) ; chop($el) ; chop($s) ; chop($p) ; chop($d) ; } # end else bond line # assume pi if no s and > 90% p- which could be dangerous!! $pi = 1 if $s == 0.0 && $p > 90.0 ; unless ( $pi ) { if ( $p > 75.0 ) { $pc = 3 ; $sc = 3.0*$s/$p ; $dc = 3.0*$d/$p ; } # end if p3 hybrid else { if ( $s != 0 ) { $sc = 1 ; $pc = $p/$s ; $dc = $d/$s ; if ( $pc > 3.0 ) { $pc = 3 ; $sc = 3.0*$s/$p ; $dc = 3.0*$d/$p ; } # end if p3d2 hybrid } # end if with s character else { $sc = 0 ; $pc = 2.0*$p/$d ; $dc = 2 ; # could be a bad assumption } # end else no s character } # end if spn hybrid } # end unless pi $sc = 0 if $sc < 0.05 ; $pc = 0 if $pc < 0.05 ; $dc = 0 if $dc < 0.05 ; $sc = substr($sc,0,4) ; $pc = substr($pc,0,4) ; $dc = substr($dc,0,4) ; print "__has $pol $el character in a " if $bond ; print "__made from a " unless $bond ; if ( $pi ) { print "p-pi orbital ($p% p"; print " $d% d" if $d>0.05 ; print ")
\n" ; } # end if pi else { if ( $sc ) { print "s$sc " unless $sc==1 ; print "s" if $sc == 1 ; } # end if s hybrid print "p$pc " if $pc ; print "d$dc " if $dc ; $hybrid = 1 ; if ( $sc ==0 && $pc ==0 ) { $hybrid = 0 } ; if ( $pc ==0 && $dc ==0 ) { $hybrid = 0 } ; print "hybrid
\n" if $hybrid ; print " orbital
\n" unless $hybrid ; } # end if not pi } # end while hybrid # now give core pairs print "

\n" ; print " -With core pairs on:" ; for ($j=1;$j<=$icr;$j++) { print "$crel[$j] " } print "-\n

\n" ; } # end for up and down print "Top of page.\n" ; # if ( $iel > 2 ) { print "

Donor Acceptor Interactions in the Best Lewis Structure

\n" ; print "The localized orbitals in your best Lewis structure\n" ; print " can interact strongly. A filled bonding or lone pair orbital can\n" ; print " act as a donor and an empty or filled bonding, antibonding, or\n" ; print " lone pair orbital can act as an acceptor. These\n" ; print " interactions can strengthen and weaken bonds. For example, a\n" ; print " lone pair donor->antibonding acceptor orbital interaction\n"; print " will weaken the bond\n" ; print " associated with the antibonding orbital. Conversly, an interaction\n" ; print " with a bonding pair as the acceptor will strengthen the bond.\n" ; print " Strong electron delocalization in your best Lewis structure will\n" ; print " also show up as donor-acceptor interactions.
\n" ; print " Interactions greater than 20 kJ/mol for bonding and lone pair\n" ; print " orbitals are listed below.\n

\n" ; while ( ) { last if /Donor/ } $_ = ; $_ = ; # skip two header lines while ( ) { if ( /within/ || /from/ ) { # pick up unit headers $x = $_ ; $_ = ; next if /None/ ; # if you want the unit headers printed out, delete the comment # # print $x ; } # end if within or from next if /^\n/ ; last if /Summary/ ; next if /CR/ || /RY/ ; # use the next line instead if you want to see core-> interactions # next if /RY/ ; ($type1,$n1,$label1,$type2,$n2,$x) = split('[()/]') ; ($j,$type1) = split('[.]',$type1) ; ($k,$type2) = split('[.]',$type2) ; ($label2,$e,$i) = split(/ {2,}/,$x) ; $e = substr($e*4.184,0,4) ; $n1 = ' the second' if $n1 == 2 ; $n1 = ' the third' if $n1 == 3 ; $n1 = $n1.'th' if $n1 > 3 ; $n1 = '' if $n1 == 1 ; $n2 = ' second' if $n2 == 2 ; $n2 = ' third' if $n2 == 3 ; $n2 = $n2.'th' if $n2 > 3 ; $n2 = '' if $n2 == 1 ; $type1 = 'antibonding' if $type1 =~ /BD\*/ ; $type1 = 'bonding' if $type1 =~ /BD/ ; $type1 = 'lone pair' if $type1 =~ /LP/ ; $type1 = 'core' if $type1 =~ /CR/ ; $type2 = 'antibonding' if $type2 =~ /BD\*/ ; $type2 = 'bonding' if $type2 =~ /BD/ ; $type2 = 'lone pair' if $type2 =~ /LP/ ; $type2 = 'core' if $type2 =~ /CR/ ; $j =~ s/ *//g ; $k =~ s/ *//g ; $label1 =~ s/ *//g ; $label2 =~ s/ *//g ; # print "The interaction of$n1 $type1 donor orbital, $j, for $label1 with \n" ; print "the$n2 $type2 acceptor orbital, $k, for $label2 is $e kJ/mol. \n" ; print "

\n" ; } # end while donor acceptor print "Top of page.\n" ; } # end if donor acceptor output # print "\n" ; print "

Molecular Orbital Energies

\n" ; print " The orbital energies are given in eV, where 1 eV=96.49 kJ/mol.\n" ; print " Orbitals with very low energy are core 1s orbitals.\n" ; print " More antibonding orbitals than you might expect are sometimes\n" ; print " listed, because d orbitals are always included for heavy\n" ; print " atoms and p orbitals are included for H atoms.\n" ; print " Up spins are shown with a ^ and down spins are shown as v.\n" ; print "Only the spin up electron orbital energies are given.\n" if $mult > 1 ; print "

\n" ; for($j=$homo+4;$j>0;$j--) { $u = '-' ; $d = '-' ; $u = '^' if $occa[$j] == 1 ; $d = 'v' if $occb[$j] == 1 ; $e = substr($moe[$j],0,5) if $moe[$j] > 0 ; $e = substr($moe[$j],0,6) if $moe[$j] < 0 ; print " $j -$u-$d- $e" ; # rough energy scaling if ( $j > 1 ) { $x=$moe[$j] ; $x=1.0 if $x<0.05 ; $x = ($moe[$j]-$moe[$j-1])/$x*100 ; $x = -$x if $x < 0 ; print "
\n" if $x > 1.0 ; print "
\n" if $x > 10.0 ; print "
\n" if $x > 100.0 ; } # else not last orbital else { print "\n" ; } } # end for all orbs # print "

\n" ; print "Top of page.\n" ; # total energy print "\n" ; print "

Total Electronic Energy

\n" ; print "The total electronic energy is a very large number, so by convention\n"; print " the units are given in atomic units, that is Hartrees (H). One\n"; print " Hartree is 2625.5 kJ/mol. The energy reference is for totally\n"; print " dissociated atoms. In other words, the reference state is a gas\n"; print " consisting of nuclei and electrons all at infinite distance from\n"; print " each other. The electronic energy includes all electric\n"; print " interactions and the kinetic energy of the electrons. This energy\n"; print " does not include translation, rotation, or vibration of the\n"; print " the molecule.\n

\n"; print " Total electronic energy = $E Hartrees\n"; # print "

\n" ; print "Top of page.\n" ; # print "

" ; print "\n"; print "-> Return to Molecular Structure Page. \n" ; print "\n"; print "-> Return to Chemistry Home Page\n" ; print "\n";