#!/usr/bin/perl
# convert mopac archive file to DeFT input file
# command line arguments are filename for archive and optimization type
# $opt='s' single point; $opt='o' optimize; $opt='r' optimize just bond lengths
# Optional arguments are $rex and $lev for use only when having SCF
# convergence problems. The bond lengths will be expanded by 1+$rex/100
# and the initial hessian guess will be dropped if $rex non-zero.
# A levelshift of $lev will be added, the SCF convergence criterion
# loosened, and the SCF mixing parameter will be changed if $lev is non-zero.
$narg = $#ARGV ;
$file = shift(@ARGV).".arc"; 	# Extract root file name into pattern
$opt = shift(@ARGV); 	# Extract runtype  into pattern
	$rex = 0 ; $lev = 0 ;
   $rex = shift(@ARGV) if $narg>1  ; # Extract %r expansion argument
   $lev = shift(@ARGV) if $narg>2 ; # Extract levelshift in eV
open(ARC,$file) || die "Can't find file $file\n";
while (<ARC>) {
	if (/FINAL/) {
		goto data ;
			}
	}
data:
$_ = scalar(<ARC>) ;
	$mult = '' ;
	$mult = 'multiplicity 2' if /DOUBLET/  ;
	$mult = 'multiplicity 3' if /TRIPLET/  ;
	$mult = 'multiplicity 4' if /QUARTET/  ;
	$charge = '' ;
	if (/CHARGE=/) {
		$charge = $_ ;
		$charge =~ s/.*CHARGE=([^ ]*).*\n/$1/ ;
		}
print $_ = scalar(<ARC>) ;
$_ = scalar(<ARC>) ;	# skip empty line
print "zmatrix angstroms\n" ;
$i = 0 ;
while (<ARC>) {            	
($sp,$el,$r,$or,$t,$oa,$p,$od,$c1,$c2,$c3,$q) = split ;
if ( $el eq '' ) {
	goto fin ;
		}
# increase bond lenghts by 10% to aid dft convergence
if ( $rex ) {
$r = substr($r*(1+$rex/100),0,5) if $r < 2.5 ;
	    }   # end if r expansion
print " $el  $c1 $r  $c2  $t  $c3  $p\n" ; 
$i++ ;
$element[$i] = $el ;
   }
fin:
print "end\n" ;
if ( $opt ne '') {
	print "runtype optimize\n" ;
	print "initial-hessian yes\n" unless $rex ;
	print "maxgeometries 30\n" ;
	print "hessian 3\n" if $opt eq 'r' ;
		}
   print "levelshift $lev\n" if $lev ;
#   loosen scf convergence criterion and change SCF mixing
   print "mixing 0.30\n" if $lev ;
   print "econverge 0.00002\n" if $lev ;
print "potential mixbecke\n" ;
print "maxiterations 51\n" ;
print "$mult\n" if $mult ne '' ;
print "charge $charge\n" if $charge ne '' ;
if ( $opt eq 'r' ){		# freeze the angles if opt=r
	print "freeze 3";
	$k = 3 ;
	for ( $j=4; $j<=$i; $j++ ) {
		$k += 2 ;
		print " ",$k ;
		$k++ ;
		print " ",$k ;
		}
			print "\n" ;
	}
print "end\n" ;

%obase1 = (
'H' , 'o-hydrogen (41)',
'Li', 'o-lithium (721/1*)',
'Be', 'o-beryllium (721/1*)',
'B' , 'o-boron (721/51/1*)',
'C' , 'o-carbon (721/51/1*)',
'N' , 'o-nitrogen (721/51/1*)',
'O' , 'o-oxygen (721/51/1*)',
'F' , 'o-fluorine (721/51/1*)',
'Al', 'o-aluminum (73111/6111/1*)',
'Si', 'o-silicon (73111/6111/1*)',
'P' , 'o-phosphorus (73111/6111/1*)',
'S' , 'o-sulfur (73111/6111/1*)',
'Cl', 'o-chlorine (73111/6111/1*)'
);
%abase1 = (
'H' , 'a-hydrogen (4;4)',
'Li', 'a-lithium (4,3;4,3)',
'Be', 'a-beryllium (5,2;5,2)',
'B' , 'a-boron (4,3;4,3)',
'C' , 'a-carbon (4,3;4,3)',
'N' , 'a-nitrogen (4,3;4,3)',
'O' , 'a-oxygen (4,3;4,3)',
'F' , 'a-fluorine (4,3;4,3)',
'AL', 'a-aluminum (5,4;5,4)',
'Si', 'a-silicon (5,4;5,4)',
'P' , 'a-phosphorus (5,4;5,4)',
'S' , 'a-sulfur (5,4;5,4)',
'Cl', 'a-chlorine (5,4;5,4)'
);
%obase2 = (
'H' , 'o-hydrogen (41/1*)',
'Li', 'o-lithium (721/1*)',
'Be', 'o-beryllium (721/1*)',
'B' , 'o-boron (7111/411/1*)',
'C' , 'o-carbon (7111/411/1*)',
'N' , 'o-nitrogen (7111/411/1*)',
'O' , 'o-oxygen (7111/411/1*)',
'F' , 'o-fluorine (7111/411/1*)',
'Al', 'o-aluminum (73111/6111/1*)',
'Si', 'o-silicon (73111/6111/1*)',
'P' , 'o-phosphorus (73111/6111/1*)',
'S' , 'o-sulfur (73111/6111/1*)',
'Cl', 'o-chlorine (73111/6111/1*)'
);
%abase2 = (
'H' , 'a-hydrogen (3,1;3,1)',
'Li', 'a-lithium (4,3;4,3)',
'Be', 'a-beryllium (5,2;5,2)',
'B' , 'a-boron (4,4;4,4)',
'C' , 'a-carbon (4,4;4,4)',
'N' , 'a-nitrogen (4,4;4,4)',
'O' , 'a-oxygen (4,4;4,4)',
'F' , 'a-fluorine (4,4;4,4)',
'AL', 'a-aluminum (5,4;5,4)',
'Si', 'a-silicon (5,4;5,4)',
'P' , 'a-phosphorus (5,4;5,4)',
'S' , 'a-sulfur (5,4;5,4)',
'Cl', 'a-chlorine (5,4;5,4)'
);

for ($j = 1; $j<=$i; $j++ ) {
$base = $abase2 { $element[$j] } ;
print "$base\n" ;
	}
for ($j = 1; $j<=$i; $j++ ) {
$base = $obase2 { $element[$j] } ;
print "$base\n" ;
	}