#!/usr/bin/perl push(@INC,"/data/htdocs/chemistry/webmo") ; require("cgi-lib.pl") ; # # make a mopac input file # command line arguments are charge, multiplicity and optimization type # $opt='s' single point; $opt='o' optimize; $opt='r' optimize just bond lengths # $opt='h' single point but calculate hessian--useful for complexes that # are unbound for AM1 but bound for DFT # followed by up to 8 elements with the connectivity implied by the order # 3 7 # | | # 4-1-2-8 # | | # 5 6 &ReadParse ; $charge = $in{'c'} ; # Extract first argument into variable $mult = $in{'m'} ; # Extract second argument into variable $opt = 'o' ; # Option that can be added later #$opt = $in{'o'} ; $el[1] = $in{'e1'} ; $el[2] = $in{'e2'} ; $el[3] = $in{'e3'} ; $el[4] = $in{'e4'} ; $el[5] = $in{'e5'} ; $el[6] = $in{'e6'} ; $el[7] = $in{'e7'} ; $el[8] = $in{'e8'} ; $file = "web".$$.".dat" ; # header for html output errors $h="Content-type: text/html\n\n

Sorry-there was an Error

"; open(DAT, "> /data/htdocs/chemistry/webmo/$file") || die "Can't redirect stdout" ; if ( $charge < -3 || $charge > 3 ) { print "$h charge out of allowed range\n" ; die ; } if ( $mult < 1 || $mult > 4 ) { print "$h multiplicity out of allowed range\n" ; die ; } if ( $opt ne 's' && $opt ne 'o' && $opt ne 'r' && $opt ne 'h' ) { print "$h optimization type not specified\n"; die ; } # load the elements $H = 0 ; $C = 0 ; $N = 0 ; $O = 0 ; $AM1 = 0 ; $MNDO = 0 ; $numsym = 31 ; # number of elements allowed + 1 for '' # this list gives the allowed elements; always end with '' @sym = ( 'H','He','Li','Be','B','C','N','O','F','Ne','Na','Mg','Al','Si', 'P','S','Cl','Ar','K','Ca','Sc','Ti','V','Cr','Mn', 'Fe','Co','Ni','Cu','Zn' ) ; # allow entry 1 if allowed 0 if not # H HeLiBeB C N O F NeNaMgAlSiP S ClArK CaScTiV CrMnFeCoNiCuZn'' @allow = (1,0,1,1,1,1,1,1,1,0,0,0,1,1,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,1) ; # atomic covalent (rep) and ionic(metals) radii T. Moeller, "Inorganic # Chemistry: A Modern Introduction," Wiley, NY, 1982 pp70 141. # H He Li Be B C N O F Ne Na Mg Al @rad = (0.37,0.5,0.68,1.1,0.8,0.77,0.74,0.74,0.71,0.65,0.97,0.66,0.51, # Si P S Cl Ar K Ca Sc Ti V Cr Mn 1.18,1.10,1.03,0.99,0.95,1.33,0.99,0.81,0.80,0.88,0.83,0.80, # Fe Co Ni Cu Zn '' 0.74,0.72,0.69,0.96,0.74,9.99 ) ; $isym = 8 ; $iatms = 0 ; for ( $j=1; $j<=$isym; $j++) { $chk = $el[$j] ; for ( $k=0; $k<=$numsym; $k++) { if ( $chk eq $sym[$k] ) { goto gotsym } } print "$h incorrect element symbol: $chk\n"; die ; gotsym: if ( ! $allow[$k] ) { print "$h element not available: $chk\n"; die } $iatms++ if $chk ne '' ; $Z[$j] = $k ; # now look for all of H C N and O $H = 1 if $chk eq 'H' ; $C = 1 if $chk eq 'C' ; $N = 1 if $chk eq 'N' ; $O = 1 if $chk eq 'O' ; $AM1 = 1 if $chk eq 'B' ; $MNDO = 1 if $chk eq 'Be' ; $MNDO = 1 if $chk eq 'Li' ; } # now rearrange atoms to make sure atoms 1 and 2 are filled # necessary for wild input for ($j=6;$j<=$isym;$j++) { if ( $el[$j] ne '' && $el[2] eq '' ) { $el[2] = $el[$j] ; $Z[2] = $Z[$j] ; $el[$j] = '' ; $Z[$j] = '' ; } } # end for terminal for ($j=3;$j<=$isym;$j++) { if ( $el[$j] ne '' && $el[1] eq '' ) { $el[1] = $el[$j] ; $Z[1] = $Z[$j] ; $el[$j] = '' ; $Z[$j] = '' ; } } # end for terminal if ( $iatms > 2 ) { unless ( $el[1] && ($el[3] || $el[4] || $el[5]) ) { # atom 1 centered $el[4] = $el[1]; $Z[4] = $Z[1] ; $el[3] = $el[7]; $Z[3] = $Z[7] ; $el[7] = '' ; $el[5] = $el[6]; $Z[5] = $Z[6] ; $el[6] = '' ; $el[1] = $el[2]; $Z[1] = $Z[2] ; $el[2] = '' ; $el[2] = $el[8]; $Z[2] = $Z[8] ; $el[8] = '' ; } # unless } # endif not diatomic # now start mopac input file print DAT " PM3" unless $AM1 || $MNDO ; print DAT " AM1" if $AM1 && !$MNDO ; print DAT " MNDO" if $MNDO ; print DAT " 1SCF" if $opt eq 's' || $opt eq 'h' ; print DAT " CHARGE=$charge" if $charge != 0 ; print DAT " DOUBLET" if $mult == 2 ; print DAT " TRIPLET" if $mult == 3 ; print DAT " QUARTET" if $mult == 4 ; print DAT " GEO-OK" ; print DAT " FORCE ISOTOPE" if $opt ne 's' ; print DAT " NOMM" if $H && $C && $N && $O ; # no mol.mechanics peptide print DAT "\n" ; print DAT " mopac input file from www\n" ; print DAT "\n" ; # Make sure we still have something in @el. Otherwise while() will run away $eel = join('',@el); if (!( $eel =~ m/\w/ )) {print "$h $eel is an empty or invalid entry\n"; die} # set up optimization flags $ropt = $aopt = $dopt = '0' ; if ( $opt eq 'o' ) { $ropt = $aopt = $dopt = '1' } if ( $opt eq 'r' ) { $ropt = '1' } # print first atom $j=1 ; while ( $el[$j] eq '' ) { $j++ } $e1 = $Z[$j] ; print DAT " $el[$j]\n" ; if ( $iatms == 1 ) { goto empty } # print second atom $j++ ; while ( $el[$j] eq '' ) { $j++ } $e2 = $Z[$j] ; $r = $rad[$e1] + $rad[$e2] ; $r = substr($r,0,5) ; print DAT " $el[$j] $r $ropt 0.000 0 0.000 0 1 0 0\n" ; if ( $iatms == 2 ) { goto empty } $jatm = 3 ; # # complete the remainder of initial fragment # get guess at hybridization $prehybrid =0 ; # previous hybrid # get number of connected atoms for ($ictr = 1;$ictr<3;$ictr++) { $icnct = 0 ; $istr = 2 ; $iend = 6 ; if ( $ictr == 2 ) { $icnct = 1 ; $istr = 6 ; $iend = 9 } for ( $k=$istr; $k<$iend;$k++) { $icnct++ if $el[$k] ne '' ; } #next neighbor k$ # hybridization exceptions @c3sp3 = ('N','O','F','P','S','Cl'); @c2sp3 = ('O','F','S','Cl'); @c2sp2 = ('N','P'); # gethybrid: { if ($icnct == 4) { $hybrid =3; last gethybrid; } if ($icnct == 3) { $hybrid =2; foreach (@c3sp3) { $hybrid=3 if $el[$ictr] eq $_ } last gethybrid; } if ($icnct == 2) { $hybrid =1; foreach (@c2sp3) { $hybrid=3 if $el[$ictr] eq $_ } foreach (@c2sp2) { $hybrid=2 if $el[$ictr] eq $_ } last gethybrid; } } # end gethybrid # geometry for initial atom center [arrays start with 0...n] @ang = (0.0,180.0,120.0,109.5) ; # starting dihedral angles 3-6 trans for sp3-sp3 and gets acetaldhyde @di2 = (0.0,0.0,180.0,180.0) ; # sp2-sp2 sp2-sp3 @di3 = (0.0,0.0,120.0,180.0) ; # sp3-sp2 sp3-sp3 $di = 0.0 ; # default starting dihedral angle $di = $di2[$hybrid] if $prehybrid == 2 ; # sp2 $di = $di3[$hybrid] if $prehybrid == 3 ; # sp3 @diinc = (0.0,0.0,180.0,120.0) ; # atom connectivity in order of fragment type; $ic1[$ictr] @ic1 = (0,1,2) ; @ic2 = (0,2,1) ; @ic3 = (0,3,3) ; $phi = $di ; $istr = $j+1 ; $iend = 6 ; if ( $ictr == 2 ) { $istr = 6 ; $iend = 9 } for ($j=$istr; $j<$iend ; $j++) { # each terminal atom $e1 = $Z[$ictr] ; $e2 = $Z[$j] ; $r = $rad[$e1] + $rad[$e2] ; $r = substr($r,0,5) ; $t = $ang[$hybrid] ; $t .= '.0' unless $t =~ /\./ ; $p = $phi ; # keep in float format $p .= '.0' unless $p =~ /\./ ; if ( $el[$j] ) { if ( $jatm != 3 ) { print DAT " $el[$j] $r $ropt $t $aopt $p $dopt", " $ic1[$ictr] $ic2[$ictr] $ic3[$ictr]\n" ; } # end if jatm not 3 else { print DAT " $el[$j] $r $ropt $t $aopt 0.000 0", " $ic1[$ictr] $ic2[$ictr] 0\n" ; } # end else jatm=3 $jatm++ ; $phi += $diinc[$hybrid] ; $phi -= 360.0 if $phi > 180.0 ; } # end unless empty } # end for terminal on atom ictr $prehybrid = $hybrid ; } # end for ictr=1..n # print out last empty line empty: print DAT "\n" ; # close(DAT) ; system("rcp /data/htdocs/chemistry/webmo/$file iris12: "); system("rsh iris12 rungms $$"); #system("rm /data/htdocs/chemistry/webmo/$file");