#!/usr/bin/perl
# 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
open(DAT, "> web$$.dat") || die "Can't redirect stdout" ;
$charge = shift(@ARGV) ; 	# Extract first argument into pattern
$mult = shift(@ARGV); 	# Extract second argument into pattern
$opt = shift(@ARGV) ;
if ( $charge < -3 || $charge > 3 ) {
		die "charge out of allowed range\n" ;
				}
if ( $mult < 1 || $mult > 4 ) {
		die "multiplicity out of allowed range\n" ;
				}
if ( $opt ne 's' && $opt ne 'o' && $opt ne 'r' && $opt ne 'h' ) {
		die "optimization type not specified\n";
				}
# 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 = $#ARGV+1 ;
$iatms = 0 ;
for ( $j=1; $j<=$isym; $j++) {
   $chk = shift(@ARGV) ;
	for ( $k=0; $k<=$numsym; $k++) {
		if ( $chk eq $sym[$k] ) { goto gotsym }
					}
			die "incorrect element symbol: $chk\n";
gotsym:
	if ( ! $allow[$k] ) { die "element not available: $chk\n" }
	$iatms++ if $chk ne '' ;
	$el[$j] = $chk ;
	$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 " 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" ;
# 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/web$$.dat iris12: ");
#system("rsh iris12 rungms $$");