#!/usr/bin/perl
use Math::Trig;
my $in_file = "Cs2SnI6.csv";
my $out_file = "Cs2SnI6";

# Read the XDACAR file
my $bandgap = 1.26;
my $energy_width = 0.001;
my @outdata_energy_minimum;
my @outdata_charge;

my $n_outdata_energies = $bandgap/$energy_width+1;

# Charge neutrality
my $cell_volume = (11.824*10**(-8))**3;
my $temperature_defect = 473.15; #unit: K.
my $temperature_defect_start = 300; #unit: K.
my $temperature_defect_end = 700; #unit: K.
my $temperature_measure = 300; #unit: K.
my $entrapy_change = 1; #unit: kBT, where KB is the boltzmann constant.
my $mass_e = 0.558; #Effective electron mass, in unit of the free electron mass, m0.
my $mass_h = 6.562; #Effective hole mass, in unit of the free electron mass, m0.
  
# Generate EF values for outdata
$outdata_energy_minimum[0][0] = "EF/eV";
$outdata_charge[0][0] = "EF/eV";
for (my $i = 1; $i < $n_outdata_energies +1; $i ++){
  $outdata_energy_minimum[$i][0] = ($i-1)*$energy_width;
  $outdata_charge[$i][0] = ($i-1)*$energy_width;
}

# Read the input data.
my @indata;
my $n_conditions; # Number of chemical potential points.
my $n_datalines; # Number of data lines in the input file. "$n_datalines - 1" is the eaxt number of data line excluing column name.
open (IN, $in_file) or die "Error: No XDATCAR file was found";
for (my $i = 0; ; $i ++){
  chomp (my $readline = <IN>);
  my @line = split (",", $readline);
  $n_conditions = @line - 3;   # 3 includes "Defect", "Charge", and "N" columns.
  for (my $j = 0; $j < @line; $j ++){ #"@line" here indicate number of data per line.
    $indata[$i][$j] = $line[$j];
    print $indata[$i][$j];
    if ($j < @line - 1){
     print ",";
    }
  }
  $n_datalines = $i;
  print "\n";
  last if (eof(IN));  
}
close (IN);

print "\nThere are ".$n_conditions." chemical potential points.\n\n";


##########################################  
##########   Data for Ploting   ##########
##########################################

# For each chemical potential.
for (my $condition = 1; $condition <= $n_conditions; $condition ++){

  my $n_defects = 0; #Number of defects considered.
  my $n_charge_states; #Number of charge states for each defect.
  # For each defect.
  for (my $i = 1; $i <= $n_datalines; $i += $n_charge_states){
    for (my $q = 1; ; $q ++){
      if ($indata[$i+$q][0] ne $indata[$i][0]){ #Determie the same defects by defect name and count the number charge states.
        $n_charge_states = $q;
        last;
      }
    }
    $n_defects += 1;
    $outdata_energy_minimum[0][$n_defects] = $indata[$i][0]; #The column name of outdata.
    $outdata_charge[0][$n_defects] = $indata[$i][0]; #The column name of outdata_charge.
    
    print  $outdata_energy_minimum[0][$n_defects].": ".$n_charge_states." charge states\n";
    
    #Obtain the minimum formation energy from various charge states for each defect.    
    for (my $j = 1; $j <= $n_outdata_energies; $j ++){
      my $minimum_energy = $indata[$i][$condition+2] + $indata[$i][1]*$outdata_energy_minimum[$j][0]; #H(EF)=H(VBM)+q*EF.
      my $minimum_charge = $indata[$i][1];
      for (my $q =1; $q < $n_charge_states; $q ++){
        if ($indata[$i+$q][$condition+2] + $indata[$i+$q][1]*$outdata_energy_minimum[$j][0] < $minimum_energy){
          $minimum_energy = $indata[$i+$q][$condition+2] + $indata[$i+$q][1]*$outdata_energy_minimum[$j][0]; #H(EF)=H(VBM)+q*EF.
          $minimum_charge = $indata[$i+$q][1];
        }
      }
      $outdata_energy_minimum[$j][$n_defects] = $minimum_energy; #Data
      $outdata_charge[$j][$n_defects] = $minimum_charge; #Data
    }
  }
  
  my $out_file_name = $out_file."-".$condition."-"."FormationEnergyMinimum.csv";
  open (OUT, ">$out_file_name");
  for (my $i = 0; $i < $n_outdata_energies+1; $i ++){
    for (my $j =0; $j <= $n_defects; $j ++){
      print OUT $outdata_energy_minimum[$i][$j];
      if ($j < $n_defects){
        print OUT ",";
      }
    }
    print OUT "\n";
  }
  close (OUT);  
  
  my $out_file_name = $out_file."-"."ChargeState.csv";
  open (OUT, ">$out_file_name");
  for (my $i = 0; $i <= $n_outdata_energies+1; $i ++){
    for (my $j =0; $j <= $n_defects; $j ++){
      print OUT $outdata_charge[$i][$j];
      if ($j < $n_defects){
        print OUT ",";
      }
    }
    print OUT "\n";
  }
  close (OUT);

  my @data_transition;
  my $n_data_transition_maximum = 0;
  for (my $i= 0; $i < $n_defects; $i ++){
    $data_transition[0][3*$i] = "EF/eV";
    $data_transition[0][3*$i+1] = $outdata_energy_minimum[0][$i+1];
    $data_transition[0][3*$i+2] = "TL";
    my $transition_count = 0;
    for (my $j = 1; $j < $n_outdata_energies; $j ++){
      if ($outdata_charge[$j][$i+1] ne $outdata_charge[$j+1][$i+1]){
        $transition_count += 1;
        $data_transition[$transition_count][3*$i] = ($outdata_energy_minimum[$j][0] + $outdata_energy_minimum[$j+1][0])/2;
        $data_transition[$transition_count][3*$i+1] = ($outdata_energy_minimum[$j][$i+1] + $outdata_energy_minimum[$j+1][$i+1])/2;
        $data_transition[$transition_count][3*$i+2] = $outdata_charge[$j][$i+1]."/".$outdata_charge[$j+1][$i+1];
        if ($transition_count > $n_data_transition_maximum){
          $n_data_transition_maximum = $transition_count;
        }
      }
    }
  }
  
  my $out_file_name = $out_file."-".$condition."-"."TransitionLevel.csv";
  open (OUT, ">$out_file_name");
  for (my $i = 0; $i <= $n_data_transition_maximum; $i ++){
    for (my $j =0; $j < 3*$n_defects; $j ++){
      print OUT $data_transition[$i][$j];
      if ($j < 3*$n_defects - 1){
        print OUT ",";
      }
    }
    print OUT "\n";
  }
  close (OUT);
}

#########################################
#########################################  
##########  Charge Neutrality  ##########
#########################################
my $out_file_name = $out_file."-T-dependence.csv";
open (OUT, ">$out_file_name");
print "Temperature";
print OUT "Temperature";
for (my $condition = 1; $condition <= $n_conditions; $condition ++){
  print ",EF/eV(".$condition."),CarrierConcentration/cm-3(".$condition.")";
  print OUT ",EF/eV(".$condition."),CarrierConcentration/cm-3(".$condition.")";
}
print "\n";
print OUT "\n";

my @outdata_T_dependence;
for (my $temperature_defect = $temperature_defect_start; $temperature_defect <= $temperature_defect_end; $temperature_defect += 25){
  print $temperature_defect;
  print OUT $temperature_defect;
  # For each chemical potential.
  for (my $condition = 1; $condition <= $n_conditions; $condition ++){
    my @data_energy; #Formation energy data as function of EF for each defect with each charge state, !!! not the minimum value.
    my @data_defect_concentration;
    my @data_charge_concentration;  
    for (my $i = 0; $i < $n_outdata_energies+1; $i ++){#Copy the Fermi energy column.
      $data_energy[$i][0] = $outdata_energy_minimum[$i][0];
      $data_defect_concentration[$i][0] = $outdata_energy_minimum[$i][0];
      $data_charge_concentration[$i][0] = $outdata_energy_minimum[$i][0];
    }
    for (my $j =1; $j < $n_datalines; $j ++){#Copy the defect name row.
        $data_energy[0][$j] = $indata[$j][0]."(".$indata[$j][1].")";
        $data_defect_concentration[0][$j] = $indata[$j][0]."(".$indata[$j][1].")";
        $data_charge_concentration[0][$j] = $indata[$j][0]."(".$indata[$j][1].")";   
    }  
    for (my $i = 1; $i < $n_outdata_energies+1; $i ++){
      for (my $j =1; $j < $n_datalines; $j ++){
        $data_energy[$i][$j] = $indata[$j][$condition+2]+$indata[$j][1]*$data_energy[$i][0];
        $data_defect_concentration[$i][$j] = ($indata[$j][2]/$cell_volume)*exp($entrapy_change - $data_energy[$i][$j]/($temperature_defect*0.00008616));
        $data_charge_concentration[$i][$j] = ($indata[$j][2]/$cell_volume)*exp($entrapy_change - $data_energy[$i][$j]/($temperature_defect*0.00008616))*$indata[$j][1]*(-1); #Negative values for electrons and positive values for holes.
      }
    }  
    # Statistical Physics of Semiconductor.
    $data_charge_concentration[0][$n_datalines]= "TotalQ-Defects";
    $data_charge_concentration[0][$n_datalines+1]= "n-StatisticalPhysics";
    $data_charge_concentration[0][$n_datalines+2]= "p-StatisticalPhysics";
    $data_charge_concentration[0][$n_datalines+3]= "Total-StatisticalPhysics";  
    # Calculate Effetive Desensity of States.
    my $m0 = 9.11*10**(-28); #unit: g.
    my $Plank = 6.626*10**(-27); #unit: gcm2/s.
    my $Bolzmann = 1.381*10**(-16); #unit: gcm2/s2/K.
    my $Bolzmann2 = 8.616*10**(-5); #unit: eV/k.
  #  my $NC = 2*(2*pi*$mass_e*$m0*$Bolzmann*$temperature_measure/$Plank/$Plank)**(3/2); #unit: cm-3.
  #  my $NV = 2*(2*pi*$mass_h*$m0*$Bolzmann*$temperature_measure/$Plank/$Plank)**(3/2); #unit: cm-3.
     my $NC = 1.8861*10**17; #unit: cm-3.
     my $NV = 8.7946*10**19; #unit: cm-3.  
  
    for (my $i = 1; $i < $n_outdata_energies+1; $i ++){
      $data_charge_concentration[$i][$n_datalines]= 0;
      for (my $j =1; $j < $n_datalines; $j ++){     
        $data_charge_concentration[$i][$n_datalines] += $data_charge_concentration[$i][$j];
      }
      $data_charge_concentration[$i][$n_datalines+1] = $NC*exp(($data_charge_concentration[$i][0] - $bandgap)/($Bolzmann2*$temperature_measure))*(-1);
      $data_charge_concentration[$i][$n_datalines+2] = $NV*exp((0 - $data_charge_concentration[$i][0])/($Bolzmann2*$temperature_measure));
      $data_charge_concentration[$i][$n_datalines+3] = $data_charge_concentration[$i][$n_datalines+1] + $data_charge_concentration[$i][$n_datalines+2];
    } 
    #Find the equilibrium Fermi energy.
    my $Fermi_level_equlibirium;
    for (my $i = 1; $i < $n_outdata_energies+1; $i ++){
      if(($data_charge_concentration[$i][$n_datalines]-$data_charge_concentration[$i][$n_datalines+3])*($data_charge_concentration[$i+1][$n_datalines]-$data_charge_concentration[$i+1][$n_datalines+3]) <= 0){
        $outdata_T_dependence[$temperature_defect - $temperature_defect_start +1][0] = $temperature_defect;
        $outdata_T_dependence[$temperature_defect - $temperature_defect_start +1][$condition] = $data_charge_concentration[$i+1][0];
        print ",".$data_charge_concentration[$i+1][0];
        print OUT ",".$data_charge_concentration[$i+1][0];
        $outdata_T_dependence[$temperature_defect - $temperature_defect_start +1][$condition + $n_conditions] = $data_charge_concentration[$i+1][$n_datalines+3];
        print ",".$data_charge_concentration[$i+1][$n_datalines+3]*(-1);
        print OUT ",".$data_charge_concentration[$i+1][$n_datalines+3]*(-1);
        last;
      }
    }
  }
  print "\n";
  print OUT "\n";
}
close (OUT);

exit;