#!/usr/bin/perl # # If you publish results please cite the following reference: # B. J. Lynch, D. G. Truhlar J. Chem. Phys. A, submitted. # # #method keywords $method = "\%mem\=200mb\n".'#AM1 scf=tight fchk'; # an @-file can be used for a user-defined basis # $ud_basis='@/home/nf/lynchb/basis/mg3s'; # or $ud_basis can be left blank for standard basis sets #$ud_basis=''; #experimental data for AE6 and BH6 benckmarks @exp_data_nospin=('6.54','19.61','10.45','12.9','3.55','17.27','322.83','192.74','102.79','705.06','633.99','1149.37'); #create g98 input files given basis and method information @input_files = &create_g98_input($method,$ud_basis); #run GAUSSIAN98 and read the output (total of 21 energies) @energies = &rung98_readg98(@input_files); #calculate atomization energies, barrier heights, AE6, and BH6 values ($ae6,$bh6,@calculated_values) = &calc_delta_energies(@energies); #print out desired information print "\nSUMMARY\n".'-------'."\n\nMETHOD\:\n"; print "$method\n\n\n"; print "\nAE6\: "; printf("%.1f",$ae6); print "\n"; print "\nBH6\: "; printf("%.1f",$bh6); print "\n"; ###END OF PROGRAM### ################# sub calc_delta_energies{ my @nrgs=@_; my @calc_vals; $calc_vals[0] = 627.5095*($nrgs[18]-$nrgs[6]-$nrgs[9] ); $calc_vals[1] = 627.5095*($nrgs[18]-$nrgs[8]-$nrgs[10] ); $calc_vals[2] = 627.5095*($nrgs[19]-$nrgs[0]-$nrgs[6] ); $calc_vals[3] = 627.5095*($nrgs[19]-$nrgs[2]-$nrgs[5] ); $calc_vals[4] = 627.5095*($nrgs[20]-$nrgs[0]-$nrgs[11] ); $calc_vals[5] = 627.5095*($nrgs[20]-$nrgs[5]-$nrgs[7] ); $calc_vals[6] = 627.5095*(4*$nrgs[0]+1*$nrgs[3] -$nrgs[12]); $calc_vals[7] = 627.5095*( 1*$nrgs[2]+1*$nrgs[3]-$nrgs[13]); $calc_vals[8] = 627.5095*( 2*$nrgs[4]-$nrgs[14]); $calc_vals[9] = 627.5095*(4*$nrgs[0]+3*$nrgs[1] -$nrgs[15]); $calc_vals[10] = 627.5095*(2*$nrgs[0]+2*$nrgs[1]+2*$nrgs[2]-$nrgs[16]); $calc_vals[11] = 627.5095*(8*$nrgs[0]+4*$nrgs[1] -$nrgs[17]); my $bar_sum=0; # for $i (0 .. 5){ # $bar_sum+=($calc_vals[$i] - $exp_data_nospin[$i])**2; # } # $bar_sum=($bar_sum/6)**(1/2); for $i (0 .. 5){ print "\n"; print abs($calc_vals[$i] - $exp_data_nospin[$i]); $bar_sum+=abs($calc_vals[$i] - $exp_data_nospin[$i]); } $bar_sum=($bar_sum/6); my $at_sum=0; # for $i (6 .. 11){ # $at_sum+=($calc_vals[$i] - $exp_data_nospin[$i])**2; # } # $at_sum=($at_sum/6)**(1/2); for $i (6 .. 11){ print "\n"; print abs($calc_vals[$i] - $exp_data_nospin[$i]); $at_sum+=abs($calc_vals[$i] - $exp_data_nospin[$i]); } $at_sum=($at_sum/6); return ($at_sum,$bar_sum,@calc_vals); } ################# sub rung98_readg98{ my @inp_files = @_; system("rm -f g98_output_files"); for $i (1 .. @inp_files){ open(FILE1,">$i.g98"); print FILE1 $inp_files[$i-1]; close FILE1; system("g98 $i.g98 $i.out"); system("cat $i.out >> g98_output_files"); open(FILE2, "){ if ($_ =~ /Total Energy/g){@temp=($_=~ /\S+/g);$pp=sprintf("%.12e",$temp[1]);$nrg[$i-1]=$pp;last;} } close FILE2; system("rm -f $i.g98"); system("rm -f $i.out"); system("rm -f Test.FChk"); } return @nrg; } ################# sub create_g98_input{ my $meth = $_[0]; my $bas = $_[1]; my @inp_files = &initialize_input; #add method keywords to beginning of the input file for $i (0 .. @inp_files-1){$inp_files[$i]=$meth."\n".$inp_files[$i]} #add user defined basis to end of input file for $i (0 .. @inp_files-1){$inp_files[$i].="\n".$bas."\n\n\n"} return @inp_files; } ################# sub initialize_input{ my @new_inputs; push (@new_inputs,' hydrogen 0 2 H '); push (@new_inputs,' carbon 0 3 C '); push (@new_inputs,' oxygen 0 3 O '); push (@new_inputs,' silicon 0 3 Si '); push (@new_inputs,' sulfur 0 3 S '); push (@new_inputs,' hydrogen molecule 0 1 H H 1 r1 r1=0.74187646 '); push (@new_inputs,' hydroxyl radical 0 2 O H 1 r1 r1=0.96889819 '); push (@new_inputs,' SH radical 0 2 S H 1 r1 r1=1.34020229 '); push (@new_inputs,' CH3 0 2 C H,1,CH H,1,CH,2,120. H,1,CH,2,120.,3,180.,0 CH=1.07731727 '); push (@new_inputs,' CH4 0 1 C H,1,RCH H,1,RCH,2,109.471221 H,1,RCH,2,109.471221,3,109.471221,1 H,1,RCH,2,109.471221,3,109.471221,-1 RCH=1.08744517 '); push (@new_inputs,' H2O 0 1 O H,1,OH H,1,OH,2,HOH OH=0.95691441 HOH=104.51706026 '); push (@new_inputs,' H2S 0 1 S 0.000000 0.000000 0.102519 H 0.000000 0.966249 -0.820154 H 0.000000 -0.966249 -0.820154 '); push (@new_inputs,' silane 0 1 Si H,1,R H,1,R,2,109.471221 H,1,R,2,109.471221,3,120.,0 H,1,R,2,109.471221,3,-120.,0 R=1.47670511 '); push (@new_inputs,' SiO 0 1 Si O,1,sio sio=1.51266699 '); push (@new_inputs,' S2 0 3 S S,1,R R=1.89259426 '); push (@new_inputs,' propyne 0 1 C 0.000000 0.000000 0.219507 C 0.000000 0.000000 1.423958 C 0.000000 0.000000 -1.243764 H 0.000000 0.000000 2.486256 H 0.000000 1.019009 -1.628154 H 0.882487 -0.509504 -1.628154 H -0.882487 -0.509504 -1.628154 '); push (@new_inputs,' glyoxal 0 1 C C,1,rcc O,1,rco,2,a1 H,1,rch,2,a2,3,180.,0 O,2,rco,1,a1,3,180.,0 H,2,rch,1,a2,3,0.,0 Variables: rcc=1.52071235 rco=1.20370135 rch=1.10078486 a1=121.36275989 a2=115.12063284 '); push (@new_inputs,' cyclobutane 0 1 C 0.000000 1.076294 0.142865 C 0.000000 -1.076294 0.142865 C -1.076294 0.000000 -0.142865 C 1.076294 0.000000 -0.142865 H 0.000000 1.979204 -0.465165 H 0.000000 1.359222 1.195822 H 0.000000 -1.979204 -0.465165 H 0.000000 -1.359222 1.195822 H -1.979204 0.000000 0.465165 H -1.359222 0.000000 -1.195822 H 1.979204 0.000000 0.465165 H 1.359222 0.000000 -1.195822 '); push (@new_inputs,' ch4oh sad pt 0 2 C -1.211487 .007968 .000407 O 1.293965 -.108694 .000133 H .009476 -.118020 .002799 H -1.525529 -.233250 1.010070 H -1.430665 1.033233 -.278082 H -1.552710 -.710114 -.737702 H 1.416636 .849894 -.000591 '); push (@new_inputs,' oh2 sad pt 0 3 H .000000 .000000 -.860287 O .000000 .000000 .329024 H .000000 .000000 -1.771905 '); push (@new_inputs,' hh2s sad pt 0 2 H 1.262097 -.220097 .000000 S .000000 .223153 .000000 H -.500576 -1.115445 .000000 H -.761521 -2.234913 .000000 '); return @new_inputs; }