#!/usr/bin/perl

#========================================
# Common Parameters
#========================================
my $Bandgap = 1.26;
my $EFWidth = 0.001;
my $NEF = $Bandgap/$EFWidth+1;
my $TM = 300;
my $EnthalpyDataFile = "Cs2SnI6.csv";
my $MaterName = "Cs2SnI6";
my $TDBEG = 300;
my $TDEND = 700;
my $TDSTEP = 25;
my $TD0 = 473.15;
my $Entrapy = 10; #unit: kBT.

#========================================
# Output Parameters
#========================================
my $OutPutEnthalpyMinimum = 1;
my $OutPutTransitionLevel = 1;
my $OutPutCHGMIN = 0;
my $OutPutDefectDensity = 1;
my $OutPutCarrierDensity = 1;
my $OutPutCarierTD = 1;

#========================================
# Read the Enthalpy Data
#========================================
my @InDataEnthalpy;
my $NChemPoint;
my $NInDataEnthalpy;
my $Flag = 0;
open (IN, $EnthalpyDataFile) or die "Error: No EnthalpyDataFile file was found";
for (my $i = 0; ; $i ++){
  chomp (my $READLINE = <IN>);
  my @InDataEnthalpyROW = split (",", $READLINE);
  $NChemPoint = @InDataEnthalpyROW - 3;
  for (my $j = 0; $j < @InDataEnthalpyROW; $j ++){ #"@InDataEnthalpyROW" here indicate number of data per InDataEnthalpyROW.
    $InDataEnthalpy[$i][$j] = $InDataEnthalpyROW[$j];
    print $InDataEnthalpy[$i][$j] if ($Flag == 1);
    print "," if ($Flag == 1 && $j < @InDataEnthalpyROW - 1);
  }
  $NInDataEnthalpy = $i;
  print "\n" if ($Flag == 1);
  last if (eof(IN));  
}
close (IN);
print "NChemPoint: ".$NChemPoint.".\n";

#========================================
# Generate EnthalpyMIN, CHGMIN and TransitionLevle
#========================================
my @OutDataEnthalpyMIN;
my @OutDataCHGMIN;
$OutDataEnthalpyMIN[0][0] = $OutDataCHGMIN[0][0] = "EF/eV";
for (my $i = 1; $i < $NEF +1; $i ++){
  $OutDataEnthalpyMIN[$i][0] = $OutDataCHGMIN[$i][0] = ($i-1)*$EFWidth;  
}
for (my $ChemPoint = 1; $ChemPoint <= $NChemPoint; $ChemPoint ++){
  my $NDefect = 0;
  my $NCHG; #Number of charge states for each defect.
  for (my $i = 1; $i <= $NInDataEnthalpy; $i += $NCHG){
    for (my $q = 1; ; $q ++){
      if ($InDataEnthalpy[$i+$q][0] ne $InDataEnthalpy[$i][0]){
        $NCHG = $q;
        last;
      }
    }
    $NDefect += 1;
    $OutDataEnthalpyMIN[0][$NDefect] = $InDataEnthalpy[$i][0]; #The column name of OutDataEnthalpyMIN.
    $OutDataCHGMIN[0][$NDefect] = $InDataEnthalpy[$i][0]; #The column name of OutDataCHGMIN.
    print  $OutDataEnthalpyMIN[0][$NDefect].": ".$NCHG." CHGs\n";
    ####: Obtain the minimum formation energy from various charge states for each defect.    
    for (my $j = 1; $j <= $NEF; $j ++){
      my $EnthalpyMIN = $InDataEnthalpy[$i][$ChemPoint+2] + $InDataEnthalpy[$i][1]*$OutDataEnthalpyMIN[$j][0]; #H(EF)=H(VBM)+q*EF.
      my $EnthalpyMINCHG = $InDataEnthalpy[$i][1];
      for (my $q =1; $q < $NCHG; $q ++){
        if ($InDataEnthalpy[$i+$q][$ChemPoint+2] + $InDataEnthalpy[$i+$q][1]*$OutDataEnthalpyMIN[$j][0] < $EnthalpyMIN){
          $EnthalpyMIN = $InDataEnthalpy[$i+$q][$ChemPoint+2] + $InDataEnthalpy[$i+$q][1]*$OutDataEnthalpyMIN[$j][0]; #H(EF)=H(VBM)+q*EF.
          $EnthalpyMINCHG = $InDataEnthalpy[$i+$q][1];
        }
      }
      $OutDataEnthalpyMIN[$j][$NDefect] = $EnthalpyMIN;
      $OutDataCHGMIN[$j][$NDefect] = $EnthalpyMINCHG;
    }
  }
  ####: Generate EnthalpyMIN File
  if ($OutPutEnthalpyMinimum == 1){
    my $OUTFile = $MaterName ."-".$ChemPoint."-"."EnthalpyMIN.csv";
    open (OUT, ">$OUTFile");
    for (my $i = 0; $i <= $NEF; $i ++){
      for (my $j =0; $j <= $NDefect; $j ++){
        print OUT $OutDataEnthalpyMIN[$i][$j];
        print OUT "," if ($j < $NDefect);
      }
      print OUT "\n";
    }
    close (OUT);  
  }
  ####: Gererate CHGMIN File
  if ($OutPutCHGMIN == 1){
    my $OUTFile = $MaterName ."-"."CHGMIN.csv";
    open (OUT, ">$OUTFile");
    for (my $i = 0; $i <= $NEF+1; $i ++){
      for (my $j =0; $j <= $NDefect; $j ++){
        print OUT $OutDataCHGMIN[$i][$j];
        if ($j < $NDefect){
          print OUT ",";
        }
      }
      print OUT "\n";
    }
    close (OUT);
  }
  ####: Gererate TransLevel File
  if ($OutPutTransitionLevel == 1){
    my @OutDataTransLevel;
    my $NOutDataTransLevelMAX = 0;
    for (my $i= 0; $i < $NDefect; $i ++){       $OutDataTransLevel[0][3*$i] = "EF/eV";
      $OutDataTransLevel[0][3*$i+1] = $OutDataEnthalpyMIN[0][$i+1];
      $OutDataTransLevel[0][3*$i+2] = "TL";
      my $TransitionCount = 0;
      for (my $j = 1; $j < $NEF; $j ++){
        if ($OutDataCHGMIN[$j][$i+1] ne $OutDataCHGMIN[$j+1][$i+1]){
          $TransitionCount += 1;
          $OutDataTransLevel[$TransitionCount][3*$i] = ($OutDataEnthalpyMIN[$j][0] + $OutDataEnthalpyMIN[$j+1][0])/2;
          $OutDataTransLevel[$TransitionCount][3*$i+1] = ($OutDataEnthalpyMIN[$j][$i+1] + $OutDataEnthalpyMIN[$j+1][$i+1])/2;
          $OutDataTransLevel[$TransitionCount][3*$i+2] = $OutDataCHGMIN[$j][$i+1]."/".$OutDataCHGMIN[$j+1][$i+1];
          if ($TransitionCount > $NOutDataTransLevelMAX){
            $NOutDataTransLevelMAX = $TransitionCount;
          }
        }
      }
    }
    my $OUTFile = $MaterName ."-".$ChemPoint."-"."TransitionLevel.csv";
    open (OUT, ">$OUTFile");
    for (my $i = 0; $i <= $NOutDataTransLevelMAX; $i ++){
      for (my $j =0; $j < 3*$NDefect; $j ++){
        print OUT $OutDataTransLevel[$i][$j];
        print OUT "," if ($j < 3*$NDefect - 1);
      }
      print OUT "\n";
    }
    close (OUT);
  }
}

#========================================
# Read POSCAR to Calculate Cell Volume
#========================================
my $CellVolume;
my @POSCAR;
my $UniversalScalingFactor;
my @Vector;
my $Flag = 0;
open (IN, "POSCAR") or die "Error: No POSCAR was found";
for (my $i = 0; ; $i ++){
  chomp ($POSCAR[$i] = <IN>);
  print $POSCAR[$i]."\n" if ($Flag == 1);
  $_ = $POSCAR[$i];
  if ($i == 1){
    if (/\s(\S+)/){
      $UniversalScalingFactor = $1;
    }
  }
  if ($i == 2){
      if (/\s(\S+)\s+(\S+)\s+(\S+)/){
        $Vector[0][0] = $1;
        $Vector[0][1] = $2;
        $Vector[0][2] = $3;
      }
  }
  if ($i == 3){
      if (/\s(\S+)\s+(\S+)\s+(\S+)/){
        $Vector[1][0] = $1;
	$Vector[1][1] = $2;
        $Vector[1][2] = $3;
      }
  }
  if ($i == 4){
      if (/\s(\S+)\s+(\S+)\s+(\S+)/){
        $Vector[2][0] = $1;
	$Vector[2][1] = $2;
        $Vector[2][2] = $3;
      }
  }  
  last if (eof(IN));
}
close (IN);
$CellVolume = 10**(-24)*$UniversalScalingFactor**3*abs(($Vector[0][1]*$Vector[1][2]-$Vector[1][1]*$Vector[0][2])*$Vector[2][0]+($Vector[0][2]*$Vector[1][0]-$Vector[1][2]*$Vector[0][0])*$Vector[2][1]+($Vector[0][0]*$Vector[1][1]-$Vector[1][0]*$Vector[0][1])*$Vector[2][2]);
print "Cell Volume: ".$CellVolume." cm-3.\n";

#========================================
# Read DOSCAR for n and p Calculations
#========================================
my @DOSCAR;
my $DOSEMAX;
my $DOSEMIN;
my $NEDOS = 10;
my $DOSWidth;
my $DOSFermi;
my $Flag = 0;
open (IN, "DOSCAR") or die "Error: No POSCAR was found";
for (my $i = 0; $i <= $NEDOS + 5; $i ++) {
  chomp ($DOSCAR[$i] = <IN>);
  print $DOSCAR[$i]."\n" if ($Flag == 1);  
  if ($i == 5){
    $_ = $DOSCAR[$i];
    if (/\s(\S+)\s+(\S+)\s+(\S+)\s+(\S+)\s+(\S+)/){
      $DOSEMAX = $1;
      $DOSEMIN = $2;
      $NEDOS = $3;
      $DOSWidth = ($1-$2)/($3-1);
      $DOSFermi = $4;
    }
  } 
}
close (IN);
print "DOSEMAX: ".$DOSEMAX."eV.\nDOSEMIN: ".$DOSEMIN."eV.\nNEDOS: ".$NEDOS.".\n"."DOSWidth: ".$DOSWidth." eV.\n";

my @TDOS;
my $Flag = 0;
for (my $i = 0; $i <= $NEDOS; $i ++){
  if ($i == 0){
    $TDOS[$i][0] = "Enegy(eV)";
    $TDOS[$i][1] = "TOS(states/cell)";
  } else{
    $_ = $DOSCAR[5+$i];
    if (/\s(\S+)\s+(\S+)\s+(\S+)/){
      $TDOS[$i][0] = $1 - $DOSFermi;
      $TDOS[$i][1] = $2;      
    }
  }
  print $TDOS[$i][0].",".$TDOS[$i][1]."\n" if ($Flag == 1);
}


#========================================
# Carrier Densities Calculated from TDOS
#========================================
my @DCarrier;
my $Flag = 1;
for (my $i = 0; $i <= $NEF; $i ++){
  if ($i == 0){
    $DCarrier[$i][0] = "EF/eV";
    $DCarrier[$i][1] = "n/cm-3";
    $DCarrier[$i][2] = "p/cm-3";
    $DCarrier[$i][3] = "Total/cm-3";
  } else{
    $DCarrier[$i][0] = ($i-1)*$EFWidth;
    my $NSUM;
    my $PSUM;
    for (my $j = 1; $j <= $NEDOS; $j ++){
      if ($TDOS[$j][0] > $Bandgap/2){
        $NSUM += $TDOS[$j][1]*$DOSWidth*(1/(1+exp(($TDOS[$j][0]-$DCarrier[$i][0])/$TM/0.00008616)))/$CellVolume*(-1);
      }
      else {
        $PSUM += $TDOS[$j][1]*$DOSWidth*(1/(1+exp(($DCarrier[$i][0]-$TDOS[$j][0])/$TM/0.00008616)))/$CellVolume;
      }
    }    
    $DCarrier[$i][1] = $NSUM;
    $DCarrier[$i][2] = $PSUM;
    $DCarrier[$i][3] = $NSUM + $PSUM;
    print $DCarrier[$i][0].",".$DCarrier[$i][1].",".$DCarrier[$i][2].",".$DCarrier[$i][3]."\n" if ($Flag == 1);    
  }
}

#========================================
# Defect and Carrier Densities
#========================================
if ($OutPutDefectDensity == 1 || $OutPutCarrierDensity == 1){
  for (my $ChemPoint = 1; $ChemPoint <= $NChemPoint; $ChemPoint ++){
    my @DataEnthalpy; #Formation energy data as function of EF for each defect with each charge state, !!! not the minimum value.
    my @DataDefectDensity;
    my @DataCarrierDensity;  
    for (my $i = 0; $i < $NEF+1; $i ++){#Copy the Fermi energy column.
      $DataEnthalpy[$i][0] = $OutDataEnthalpyMIN[$i][0];
      $DataDefectDensity[$i][0] = $OutDataEnthalpyMIN[$i][0];
      $DataCarrierDensity[$i][0] = $OutDataEnthalpyMIN[$i][0];
    }
    for (my $j =1; $j < $NInDataEnthalpy; $j ++){#Copy the defect name row.
        $DataEnthalpy[0][$j] = $InDataEnthalpy[$j][0]."(".$InDataEnthalpy[$j][1].")";
        $DataDefectDensity[0][$j] = $InDataEnthalpy[$j][0]."(".$InDataEnthalpy[$j][1].")";
        $DataCarrierDensity[0][$j] = $InDataEnthalpy[$j][0]."(".$InDataEnthalpy[$j][1].")";   
    }
    $DataCarrierDensity[0][$NInDataEnthalpy]= "Total-Q-Defects";
    for (my $i = 1; $i <= $NEF; $i ++){
      $DataCarrierDensity[$i][$NInDataEnthalpy]= 0;
      for (my $j =1; $j < $NInDataEnthalpy; $j ++){
        $DataEnthalpy[$i][$j] = $InDataEnthalpy[$j][$ChemPoint+2]+$InDataEnthalpy[$j][1]*$DataEnthalpy[$i][0];
        $DataDefectDensity[$i][$j] = ($InDataEnthalpy[$j][2]/$CellVolume)*exp($Entrapy - $DataEnthalpy[$i][$j]/($TD0*0.00008616));
        $DataCarrierDensity[$i][$j] = ($InDataEnthalpy[$j][2]/$CellVolume)*exp($Entrapy - $DataEnthalpy[$i][$j]/($TD0*0.00008616))*$InDataEnthalpy[$j][1]*(-1); #Negative values for electrons and positive values for holes.
        $DataCarrierDensity[$i][$NInDataEnthalpy] += $DataCarrierDensity[$i][$j];
      }
    }
    if ($OutPutDefectDensity == 1){
      my $OUTFile = $MaterName ."-".$ChemPoint."-"."DefectDensity.csv";
      open (OUT, ">$OUTFile");
      for (my $i = 0; $i <= $NEF; $i ++){
        for (my $j =1; $j < $NInDataEnthalpy; $j ++){
          print OUT $DataDefectDensity[$i][$j];
          print OUT "," if ($j < $NInDataEnthalpy);
        }
        print OUT "\n";
      }
      close (OUT); 
    }
    
    if ($OutPutCarrierDensity == 1){
      my $OUTFile = $MaterName ."-".$ChemPoint."-"."CarrierDensity.csv";
      open (OUT, ">$OUTFile");
      for (my $i = 0; $i <= $NEF; $i ++){
        for (my $j =1; $j <= $NInDataEnthalpy; $j ++){
          print OUT $DataCarrierDensity[$i][$j];
          print OUT "," if ($j < $NInDataEnthalpy);
        }
        print OUT "\n";
      }
      close (OUT);
    }
  }
}

#========================================
# Self-Consistency by Charge Neutrality
#========================================
if ($OutPutCarierTD == 1){
  my $OUTFile = $MaterName ."-Carrier-TD.csv";
  open (OUT, ">$OUTFile");
  print "Temperature";
  print OUT "Temperature";
  for (my $ChemPoint = 0; $ChemPoint <= $NChemPoint; $ChemPoint ++){
    if ($ChemPoint == 0){
      print "Temperature";
      print OUT "Temperature";
    } else{
      print ",EF/eV(".$ChemPoint."),CarrierConcentration/cm-3(".$ChemPoint.")";
      print OUT ",EF/eV(".$ChemPoint."),CarrierConcentration/cm-3(".$ChemPoint.")";
    }
  }
  print "\n";
  print OUT "\n";
  my @OutDataCarrierDensityTD;
  for (my $TD = $TDBEG; $TD <= $TDEND; $TD += $TDSTEP){
    print $TD;
    print OUT $TD;
    for (my $ChemPoint = 1; $ChemPoint <= $NChemPoint; $ChemPoint ++){
      my @DataEnthalpy; #Formation energy data as function of EF for each defect with each charge state, !!! not the minimum value.
      my @DataDefectDensity;
      my @DataCarrierDensity;  
      for (my $i = 0; $i < $NEF+1; $i ++){#Copy the Fermi energy column.
        $DataEnthalpy[$i][0] = $OutDataEnthalpyMIN[$i][0];
        $DataDefectDensity[$i][0] = $OutDataEnthalpyMIN[$i][0];
        $DataCarrierDensity[$i][0] = $OutDataEnthalpyMIN[$i][0];
      }
      for (my $j =1; $j < $NInDataEnthalpy; $j ++){#Copy the defect name row.
          $DataEnthalpy[0][$j] = $InDataEnthalpy[$j][0]."(".$InDataEnthalpy[$j][1].")";
          $DataDefectDensity[0][$j] = $InDataEnthalpy[$j][0]."(".$InDataEnthalpy[$j][1].")";
          $DataCarrierDensity[0][$j] = $InDataEnthalpy[$j][0]."(".$InDataEnthalpy[$j][1].")";   
      }
      $DataCarrierDensity[0][$NInDataEnthalpy]= "Total-Q-Defects";
      for (my $i = 1; $i <= $NEF; $i ++){
        $DataCarrierDensity[$i][$NInDataEnthalpy]= 0;
        for (my $j =1; $j < $NInDataEnthalpy; $j ++){
          $DataEnthalpy[$i][$j] = $InDataEnthalpy[$j][$ChemPoint+2]+$InDataEnthalpy[$j][1]*$DataEnthalpy[$i][0];
          $DataDefectDensity[$i][$j] = ($InDataEnthalpy[$j][2]/$CellVolume)*exp($Entrapy - $DataEnthalpy[$i][$j]/($TD*0.00008616));
          $DataCarrierDensity[$i][$j] = ($InDataEnthalpy[$j][2]/$CellVolume)*exp($Entrapy - $DataEnthalpy[$i][$j]/($TD*0.00008616))*$InDataEnthalpy[$j][1]*(-1); #Negative values for electrons and positive values for holes.
          $DataCarrierDensity[$i][$NInDataEnthalpy] += $DataCarrierDensity[$i][$j];
        }
      }
      ####: Self-Consistency
      for (my $i = 1; $i <= $NEF; $i ++){
        if(($DataCarrierDensity[$i][$NInDataEnthalpy]-$DCarrier[$i][3])*($DataCarrierDensity[$i+1][$NInDataEnthalpy]-$DCarrier[$i+1][3]) <= 0){
          $OutDataCarrierDensityTD[$TD - $TDBEG +1][0] = $TD;
          $OutDataCarrierDensityTD[$TD - $TDBEG +1][$ChemPoint] = $DataCarrierDensity[$i+1][0];
          print ",".$DataCarrierDensity[$i+1][0];
          print OUT ",".$DataCarrierDensity[$i+1][0];
          $OutDataCarrierDensityTD[$TD - $TDBEG +1][$ChemPoint + $NChemPoint] = $DataCarrierDensity[$i+1][$NInDataEnthalpy+3];
          print ",".$DCarrier[$i+1][3]*(-1);
          print OUT ",".$DCarrier[$i+1][3]*(-1);
          last;
        }
      }
    }
    print "\n";
    print OUT "\n";
  }
  close (OUT);
}

exit;