package oniom_ct;
use strict;
use warnings;
use base 'Exporter';
our @EXPORT=qw(oniom_ct);

sub oniom_ct{
    use Raghavachari_Ext_Pack_utilities;
    use strict;
    use warnings;
    my @options=@_;
    my $spin = pop(@options);
    my $charge = pop(@options);
    my $DoDeriv = pop(@options);
    my $NAtoms = pop(@options);
    my $in_geom = pop(@options);
    my $message = pop(@options);
    my $output = pop(@options);
    my $input = pop(@options);
    my $layer = pop(@options);
    

    #   
    #   The outline is as follows:
    #       1) Submit the RL sub-calculation.
    #       2) Submit the first ML sub-calculation with standard link-atom nuclear charges.
    #       3) Read ML output, obtain regional charges, and increment the link-atom nuclear charges.
    #       4) Submit the second ML sub-calculation, repeat step 5), and then check for convergence, 
    #           for example, do the ML regional charges match the RL regional charges?
    #       5) Iterate till converged.
    #       7) Submit the MH sub-calculation with the optimized link-atom nuclear charges.
    #       8) Integrate the ONIOM-CT energy.
    #
    
    
    #
    #   Parse input keywords and set default values.
    my %opt=();
    foreach (@options){
        if($_=~m/^-(.*)==(.*)$/i){
            $opt{"$1"}="$2";
        }else{
            die" Unknown option: $_\n";
        };
    };
    my $initialstep=0;
    unless($opt{initialstep}){
        $initialstep=0.0015;
    }else{
        $initialstep=$opt{"initialstep"};
    };
    
    #   
    #   Do Forces?
    $opt{method}="$opt{method} units=au";
    if($opt{grad}){
        if($opt{grad}==1){
            $opt{method}="$opt{method} force";
        }elsif($opt{grad}==2){
            die "Can't do force constants yet!\n";
        }else{die "Bad value for option 'grad'=$opt{grad}}/n"
        };
    }else{
        $opt{grad}=0;
    };
    
    print "====================================================================\n"; 
    print " In ONIOM-CT module: Geometry picked up from Gau-.EIn file:\n";
    print "         Current implementation requires the -chk== option to \n";
    print "         have the same name as the root filename.\n";
    print "====================================================================\n"; 
    open(MESSAGE,">>","$message") or die "Couldn't open MESSAGE\n";


    #   Define path for gt
    if($opt{execute}=~m/^gt$/){
        $opt{execute} = "gdv -exedir=/home/nmayhall/gauss-h08:/home/nmayhall/gauss-h08/exe-dir:/gaussian/gdv-h08/gdv/bsd:/gaussian/gdv-h08/gdv/local:/gaussian/gdv-h08/gdv/extras:/gaussian/gdv-h08/gdv";
    };
    
    #
    #   What population analysis to use?
    my $pop="";
    unless($opt{population}){
        $pop="mull";
    }else{
        $pop=$opt{population};
    };
    if($pop=~m/^lowd/i){
        $opt{method}="$opt{method} iop(6/80=1)";
    };

    #
    #   Set up RWF file to obtain info about geom step number and other external info.
    open(RWF,">>","$opt{chk}_rwf.txt") or die "Couldn't open RWF\n";
    my $geom_iter=0;
    while(<RWF>){
        if($_ =~ m/^geometry_iteration_number=(\d+)$/){
            $geom_iter=$1;
        };
    };
    $geom_iter=$geom_iter+1;
    print RWF "geometry_iteration_number=$geom_iter\n";

    #
    #   Print out option list to log file.
    foreach (sort(keys %opt)){
        printf " %20s = %-25s\n",$_,$opt{$_};
    };
    
    
    my @in_geom=split(";",$in_geom);
    #   Ensure that the input geometry has the correct dimensions.
    unless(($#in_geom+1) % 4 == 0){die " Problem with input geometry:$#in_geom\n"};
    my $i = 0;
    

    #
    my $geom=();
    my $route=();
    my $title=();
    my $add_info=();
    #  Jovan's options
    if(not $opt{zmat}){
        $opt{zmat}="0";
    };
    if(not $opt{lacopt}){
        $opt{lacopt}="0";
    };        
    if($opt{zmat}==1){
      ($route,$title,$geom,$add_info)=parse_inputfile_zmat("$NAtoms", "$opt{chk}",";","$in_geom");
    }elsif($opt{zmat}==0){
      ($route,$title,$geom,$add_info)=parse_inputfile("$opt{chk}",";");
    }else{
      print "Warning Error: Check the zmat criterion\n";	    
    }
    #
    my ($oniom_layer)=parse_geom($geom,";");
    my @oniom_layer=split(";",$oniom_layer);
    @oniom_layer=("empty",@oniom_layer);
    print "WHAT IS LAYER @oniom_layer\n";
    unless($#oniom_layer==$NAtoms){die "Wrong number of elements in \@oniom_layer\n"};

    my @linkatom_info=();
    my @real_info=();
    my @model_info=();
    #
    #   Fill @info_* arrays 
    #       These arrays are used throughout the program by returning certain 
    #       information about the input arguement.
    #
    #       i.e.
    #       for atom center $i
    #          
    #           $info_linkatom_back[$i]
    #
    #       returns the atomic number of the link-atom support,
    #       where $i is the link atom host. 
    #        
    #
    my @info_id =();
    my @info_linkatom =();
    my @info_la_h2s =();
    my @info_linkatom_back =();
    my @info_linkatom_capID =();
    my @info_linkatom_scale =();
    my @info_model =();
    my @info_atom_ssystems =();
    #   Since these arrays correspond to center number, all zero elements are set to "empty".
    $info_id[0]="empty";
    $info_la_h2s[0]="empty";
    $info_linkatom[0]="empty";
    $info_linkatom_back[0]="empty";
    $info_linkatom_capID[0]="empty";
    $info_linkatom_scale[0]="empty";
    for($i=1;$i<=($NAtoms);$i++){
        my $j=($i-1)*4;
        $info_id[$i]=$in_geom[$j];
        if($oniom_layer[$i]=~m/^Low/){
            $info_model[$i]=1            
        }else{
            $info_model[$i]=0            
        }
        $info_la_h2s[$i]=0;
        $info_linkatom[$i]=0;
        $info_linkatom_back[$i]=0;
        $info_linkatom_capID[$i]=0;
        $info_linkatom_scale[$i]=0;
        $info_atom_ssystems[$i]="r";
    };
    for($i=1;$i<=($NAtoms);$i++){
        if($info_model[$i]==0){
            $info_atom_ssystems[$i]="$info_atom_ssystems[$i]".";m";
        }
    };
    #   
    #   Fill link atom @info arrays.
    #
    for($i=1;$i<=($NAtoms);$i++){
        if($oniom_layer[$i]=~m/^Low_with_(\S+)_link_to_(\d+)_by_(\d+\.\d+)$/){
            $info_linkatom[$2]=$i;
            $info_linkatom_back[$i]=$2;
            $info_la_h2s[$i]="$info_la_h2s[$i]".";$2";
            $info_linkatom_capID[$2]=$1;
            $info_linkatom_scale[$i]=$3;
        };
    };
    print "--Link-Atom Info:\n";
    for($i=1;$i<=($NAtoms);$i++){
        if($oniom_layer[$i]=~m/^Low_with_(\S+)_link_to_(\d+)_by_(\d+\.\d+)$/){
            my @temp = split(";",$info_la_h2s[$i]);
            shift(@temp);
            foreach my $s (@temp){
                my @temp1 = ($i,$info_id[$i],$s,$info_id[$s]);
                my @temp2 = ($info_linkatom_capID[$info_linkatom_back[$i]],$info_linkatom_scale[$i]);
                printf "Atom %5i(ID=%3i) is a link-atom host connected to %5i(ID=%3i) replaced by a %5s link-atom: scaled by %10f\n",@temp1,@temp2;
            };
            
        };
    };
    #
    #   Fill geometry arrays.
    #
    #   For link atoms, use the following definitions of terms:
    #
    #       |--C1-C2--|   -->   |--C1-H
    #       
    #       C1 == link-atom support.
    #       C2 == link-atom host.
    #       H  == link-atom.
    #

    #   Initiate geom_real and geom_model hashes.
    my (@geoms_real,@geoms_model)=();
    $i = 0;
    foreach my $k (@model_info){
        $i =$i+1;
        $geoms_model[$i]=";";
    };
    
    #   Initiate %geoms hash.
    my %geoms=();
    my @sys_unique = qw(rl ml mh);
    foreach my $sys (@sys_unique){
        my $method = "";
        if($sys=~m/^rl$/){
            $method = $opt{low};
        }elsif($sys=~m/^ml$/){
            $method = "$opt{low} massage";
        }elsif($sys=~m/^mh$/){
            $method = "$opt{high} massage";
        }else{die "Confused \$sys=$sys\n";}
        open (HANDLE,">","$opt{chk}"."_$sys.gjf") or die "Couldn't open file:$opt{chk}"."_$sys.gjf\n";
        print HANDLE "%chk=$opt{chk}"."_$sys.chk\n%mem=$opt{mem}\n%nproc=$opt{nproc}\n#p $method $opt{method}\nnosymm\n\ntest - NoSymm still required.\n\n";
        if(($opt{charge_mult_rm}) and ($sys !~ m/^f$/)){
            print HANDLE "$opt{charge_mult_rm}\n";
        }else{
            print HANDLE "$charge $spin\n";
        };
        $geoms{$sys}=";";
        close(HANDLE);
    };
   
    #   1)
    #   Set up the RL, ML(0), and MH(0) calculations. 
    #
     
    #   Fill %geoms hash with geometries for each subsystem.   
    for($i=1;$i<=($NAtoms);$i++){
        my $j=($i-1)*4;
        foreach my $sys (@sys_unique){
            open (HANDLE,">>","$opt{chk}"."_$sys.gjf") or die "Couldn't open file\n";
            my $atom="";
            my $x=$in_geom[$j+1];
            my $y=$in_geom[$j+2];
            my $z=$in_geom[$j+3];
            if($sys=~m/^rl$/){
                $atom=$in_geom[$j];
            }elsif($sys=~m/^m(l|h)$/){
                if($oniom_layer[$i]=~m/^High$/){
                    $atom=$in_geom[$j];
                }elsif($oniom_layer[$i]=~m/^Low$/){
                    $atom="Bq";
                }elsif($oniom_layer[$i]=~m/^Low_with_(\S+)_link_to_(\d+)_by_(\d+\.\d+)$/){
                    $atom=$1;
                    my $host = $i;
                    my $supp = $2;
                    my $id=$info_linkatom_capID[$supp];
                    my $scale = $info_linkatom_scale[$i];
                    if($scale==0){die "Link atom scale factors must be defined in input file\n"};
                    #   Scale Link-atom Coordinates.
                    $x=$in_geom[($supp-1)*4+1];
                    $y=$in_geom[($supp-1)*4+2];
                    $z=$in_geom[($supp-1)*4+3];
                    $x=$x+$scale*($in_geom[($host-1)*4+1]-$x);
                    $y=$y+$scale*($in_geom[($host-1)*4+2]-$y);
                    $z=$z+$scale*($in_geom[($host-1)*4+3]-$z);
                };
            };
            printf HANDLE  "%8s %20.12f %20.12f %20.12f\n",$atom,$x,$y,$z;
        };
    };


    #   Add blank line to end of each input file.
    foreach my $sys (@sys_unique){
        open (HANDLE,">>","$opt{chk}"."_$sys.gjf") or die "Couldn't open file: $opt{chk}"."_$sys.gjf\n";
        print HANDLE "\n";
        close HANDLE;
    };
    
    #   Add massage input 
    for($i=1;$i<=($NAtoms);$i++){
        my $j=($i-1)*4;
        if($oniom_layer[$i]=~m/^Low_with_(\S+)_link_to_(\d+)_by_(\d+\.\d+)$/){
            foreach my $sys (qw(ml mh)){
                open (HANDLE,">>","$opt{chk}"."_$sys.gjf") or die "Couldn't open file: $opt{chk}"."_$sys.gjf\n";
                printf HANDLE "%8i%8s%20.12f\n",$i,"NUC",1.0;
                close HANDLE;
            };
        };
    };

    #   Add blank line to end of each input file.
    foreach my $sys (qw(ml mh)){
        open (HANDLE,">>","$opt{chk}"."_$sys.gjf") or die "Couldn't open file: $opt{chk}"."_$sys.gjf\n";
        print HANDLE "\n";
        close HANDLE;
    };

    
    #   2)  
    #   Submit the RL Calculation.
    #
    
    my $new="nicknewline";
    my (%energy,%nbasis,%atomicID,%carts,%TVcarts,%rmsforce,%cartgradient,%charges,%dipmom)=();
    unless($opt{nogdv}==1){
        system("$opt{execute} $opt{chk}_rl.gjf");
    };
    system("formchk $opt{chk}_rl.chk");
    ($energy{"rl"},$nbasis{"rl"},$atomicID{"rl"},$carts{"rl"},$TVcarts{"rl"},$rmsforce{"rl"},$cartgradient{"rl"},$dipmom{"rl"},$charges{"rl"},,,)=parse_fchk("$opt{chk}_rl");
    my @cartgradient    =split(",",$cartgradient{rl});
#   Alex: let's get charges from RL log file,
#   if this is Lowdin calculation.    
    if($pop=~m/^lowd/i){
        ($charges{rl}) = get_lowdin_charges("$opt{chk}_rl.log");
    };
    my @charges_rl    =split(",",$charges{rl});

    $cartgradient{"rl"}="";
    foreach(@cartgradient){
        if($_=~m/^(-?\d\.\d+)E(-?\+?\d+)/){
            $_ = $1*10**$2;
            $cartgradient{"rl"}="$cartgradient{rl}"."$_;";
        }else{die "Couldn't read number format for cartgradient.\n"};
    };
#    my ($del,$energy_rl,$charges_rl)=ct_oniom_analysis("$opt{chk}"."_rl","1",$pop,$new);
    my $del =();

    my($qRL_A)=0;
    my $count=-1;

    my $counter = 0;
    foreach(@oniom_layer){
        $counter = $counter + 1;
        if($_ =~ m/^High$/){
            my $temp_charge = $charges_rl[$count];
            $qRL_A = $qRL_A  + $temp_charge;
        };
        $count = $count+1;
    };
    
    printf MESSAGE "%-17s %10s =   %11.7f   %10s =   %14.10f\n","Sub-System RL","Energy:",$energy{rl},"Charge:",$qRL_A;
    
    #   4)
    #   Run the initial ML input file.  
    #

    unless($opt{nogdv}==1){
        system("$opt{execute} $opt{chk}_ml.gjf");
    };
    system("formchk $opt{chk}_ml.chk");
    ($energy{"ml"},$nbasis{"ml"},$atomicID{"ml"},$carts{"ml"},$TVcarts{"ml"},$rmsforce{"ml"},$cartgradient{"ml"},$dipmom{"ml"},$charges{"ml"},,,)=parse_fchk("$opt{chk}_ml");
    
    @cartgradient    =split(",",$cartgradient{ml});
    $cartgradient{"ml"}="";
    foreach(@cartgradient){
        if($_=~m/^(-?\d\.\d+)E(-?\+?\d+)/){
            $_ = $1*10**$2;
            $cartgradient{"ml"}="$cartgradient{ml}"."$_;";
        }else{die "Couldn't read number format for cartgradient.\n"};
    };
    
    if($pop=~m/^lowd/i){
        ($charges{ml}) = get_lowdin_charges("$opt{chk}_ml.log");
    };
    my @charges_ml    =split(",",$charges{ml});
    
    my $resid=1.0;
    my $N=0;
    my ($Echarge,$Nchange,@pcharge,@qML_A,$pcharge)=();
    $pcharge[0]=1.0;
    $pcharge[1]=$initialstep;
    $pcharge[1]=$pcharge[0]+$initialstep;

    #
    #   Iterate to convergence.
    my $conv = 0.00001;
    if($opt{converge_charge}){
        $conv = $opt{converge_charge};
    };
    my $qML_A=0;
    my ($energy_ml,$charges_ml)=();
    while($resid >= $conv){
        if($N gt 1){
            my $denom_test=abs($qML_A[$N-2]-$qML_A[$N-1]);
            if($denom_test ne 0){
                $pcharge[$N] = $pcharge[$N-1]+($pcharge[$N-2]-$pcharge[$N-1])*($qRL_A-$qML_A[$N-1])/($qML_A[$N-2]-$qML_A[$N-1]);
            }else{
                $pcharge[$N]=$pcharge[$N-1];
            }
        };
        #
        # Prepare the ML[$N] input file.
        #
        $pcharge[$N] = sprintf("%9.12f",$pcharge[$N]);
        open(ML,"$opt{chk}_ml.gjf") or die;
        my $temp_mlinput=join("",<ML>);
        close(ML);
        if($N gt 1){
            # 
            # Adjust charge in multiplicity if link atom picks up 
            # more than half an electron.
            #
            if($temp_mlinput=~m/\n\s*(-?\d+)(\s+|,)(\d+)\s*\n/){
              $Echarge = $1;
            };
            $Echarge = update_molecular_charge($Echarge,$pcharge[$N-1],$pcharge[$N]);

            $temp_mlinput=~s/(\n\s*)(-?\d+)(\s+|,\d+\s*\n)/$1 $Echarge $3/;
        };
        $temp_mlinput=~s/NUC\s+-?\d+\.\d+/NUC $pcharge[$N]/g;
        if($N == 1){
            $temp_mlinput=~s/massage/massage guess=read/g;
        };
        open(ML2INPUT,">","$opt{chk}_ml.gjf");
        print ML2INPUT "$temp_mlinput\n";
        # 
        # Run gdv.
        #
        print "   **Submitting ML job:    $opt{chk}_ml.gjf:    Iteration $N**\n"; 
        unless($opt{nogdv}){
            system("$opt{execute} $opt{chk}_ml.gjf");
        };
        system("formchk $opt{chk}_ml.chk");
        
        # 
        # Read ML[$N] output. Get energy and model system charge (qML_A).
        #
    
        ($energy{"ml"},$nbasis{"ml"},$del,$del,$del,$del,$del,$del,$charges{"ml"},,,)=parse_fchk("$opt{chk}_ml");
        if($pop=~m/^lowd/i){
            ($charges{ml}) = get_lowdin_charges("$opt{chk}_ml.log");
        };  
        my @charges_ml    =split(",",$charges{ml});
        
        $qML_A=0;
        $count=-1;
        foreach(@oniom_layer){
            if($_ =~ m/^High$/){
                my $temp_charge = $charges_ml[$count];
                $qML_A = $qML_A  + $temp_charge;
            };
            $count = $count+1;
        };
        $qML_A[$N]=$qML_A;
        
        #
        # Check for convergence.
        #
        
        if($N gt 1){
            $resid=abs($qML_A[$N]-$qML_A[$N-1]);
        }else{
            $resid=abs($qRL_A-$qML_A[$N]);
        };
        $pcharge=$pcharge[$N];
        printf MESSAGE "%-15s %12s =   %11.8f   %10s =   %12.8f   %10s =   %12.8f   %10s =   %14.10f\n","Sub-System ML.$N","Energy:",$energy{ml},"Charge:",$qML_A,"Residual:",$resid,"NCharge:",$pcharge[$N];
	if($opt{lacopt}==1){last;}
        $N=$N+1;
    };
   
    #
    #   Execute MH calculation with optimized PCharge
    open(MHINPUT,"$opt{chk}_mh.gjf");
    chomp(my @mhinput = <MHINPUT>);
    close(MHINPUT);
    my $mhinput = join($new,@mhinput);
    $mhinput=~s/$new/\n/g;
    $pcharge = sprintf("%9.12f",$pcharge);
    my $temp_mhinput=$mhinput;
    $Nchange =  $pcharge[$N-1] - $pcharge[0];

    if(abs($Nchange) >= 0.5){
        # 
        # Run a separate MH job with charge 1 on the link atom 
        # to obtain the correct guess. This is necessary because
        # links 301 and 401 behave differently with respect to
        # the high charge on the link atom. For example,
        # when that charge increases from 1.0 to 1.55, 
        # the link 301 adds an extra electron
        # to the molecule, but link 401 does not.
        #
        open(MHINPUT,">","$opt{chk}_mh.gjf");
        print MHINPUT "$temp_mhinput\n";
        close(MHINPUT);
    
        print "   **Submitting a preliminary MH job:    $opt{chk}_mh.gjf\n"; 
        unless($opt{nogdv}){
            system("gdv $opt{chk}_mh.gjf");
        };
        # 
        # Adjust charge and multiplicity if link atom picks up 
        # more than half an electron.
        #
        if($temp_mhinput=~m/\n\s*(-?\d+)(\s+|,)(\d+)\s*\n/){
            $Echarge = $1;
        };
        if($Nchange <= -0.5){
            $Echarge = $Echarge-1;
        }elsif($Nchange >= 0.5){
            $Echarge = $Echarge+1;
        };
        $temp_mhinput=~s/(\n\s*)(-?\d+)(\s+|,\d+\s*\n)/$1 $Echarge $3/;
        $temp_mhinput=~s/massage/massage guess=read/g;
    };
    $temp_mhinput=~s/NUC\s+-?\d+\.\d+/NUC $pcharge/g;
    open(MHINPUT,">","$opt{chk}_mh.gjf");
    print MHINPUT "$temp_mhinput\n";
    
    print "   **Submitting MH job:    $opt{chk}_mh.gjf\n"; 
    unless($opt{nogdv}){
        system("gdv $opt{chk}_mh.gjf");
    };
    
    # 
    # Read MH output. Get energy and model system charge (qMH_A).
    #
    
    system("formchk $opt{chk}_mh.chk");
    ($energy{"mh"},$nbasis{"mh"},$del,$del,$del,$del,$del,$del,$charges{"mh"},,,)=parse_fchk("$opt{chk}_mh");
    if($pop=~m/^lowd/i){
        ($charges{mh}) = get_lowdin_charges("$opt{chk}_mh.log");
    };
    my @charges_mh    =split(",",$charges{mh});

    my $qMH_A=0;
    $count=-1;
    foreach(@oniom_layer){
        if($_ =~ m/^High$/){
            my $temp_charge = $charges_mh[$count];
            $qMH_A = $qMH_A  + $temp_charge;
        };
        $count = $count+1;
    };
    printf MESSAGE "%-15s %12s =   %11.8f   %10s =   %12.8f   %27s   %10s =   %14.10f\n","Sub-System MH","Energy:",$energy{mh},"Charge:",$qMH_A,"","NCharge:",$pcharge;
    
    #
    #   6)
    #   Integrate ONIOM-CT and print out final data.
    #
    
    my $oniomct_energy = $energy{rl} + $energy{mh} - $energy{ml};
    print MESSAGE "\n\n";
    printf MESSAGE "%-10s =   %15.11f\n","ONIOM-CT",$oniomct_energy;
    printf MESSAGE "%-10s =    %11.9f\n","PCharge",$pcharge;

    
     
    my $E =$oniomct_energy;
    my @grad=();
    
    printf "SCF Done Energy= %.18f\n",$E;
    close(MESSAGE); 
    my $grad = join(";",@grad); 
    return($E,$grad);
};