import java.io.*;
public class DoubleBondOxygenPassivation {
  public static void main(String argv[]){

    int i, j, k, m, natom;
    int sum_neighbor;

    FileInputStream fis;
    FileOutputStream Table_out, DataOut_ordered, Jmol_out;

    PrintStream prn_table, prn_jmol; 

    String atom_name;
    double [][] all_atoms;
    int [][] table_neighbor;
    int [] table_amount;
    int num_of_hydrogen;
    int num_of_oxygen;
    int num_of_silicon;
    double lowest_distance = 1.0e+20;
    double distance, norm, scalar_product;
    double [] tempo, vector, vectorA, vectorB, vectorC;

    double SiO_bondlength = 1.65;
    double SiH_bondlength = 1.50;
    double OH_bondlength =  0.96;

    tempo   = new double [3];
    vector  = new double [3];
    vectorA = new double [3];
    vectorB = new double [3];
    vectorC = new double [3];

    try {    
      fis = new FileInputStream("quantumdot.dat");

      Table_out = new FileOutputStream("NNtable_WITH_oxygen.list");
      Jmol_out = new FileOutputStream("oxygen_passivated.xyz");
      DataOut_ordered = new FileOutputStream("WITH_oxygen.dat");

      DataInputStream dis = new DataInputStream(fis);
      prn_table = new PrintStream( Table_out );
      prn_jmol = new PrintStream( Jmol_out );

      DataOutputStream Data_out = new DataOutputStream( DataOut_ordered );
      natom = dis.readInt();
      System.out.println(natom);

      all_atoms = new double [natom][3];
      table_amount   = new int [natom];
      table_neighbor = new int [natom][4];

      for( i=0; i < natom; i++) {
          atom_name = dis.readUTF();
          all_atoms[i][0] = dis.readDouble();
          all_atoms[i][1] = dis.readDouble();
          all_atoms[i][2] = dis.readDouble();
      }

      // find lowest distance
      for( i=0; i < natom; i++) {

        tempo[0] = all_atoms[i][0];
        tempo[1] = all_atoms[i][1];
        tempo[2] = all_atoms[i][2];

        for( j=0; j < natom; j++) {
          if (j != i) {
            distance = Math.pow(Math.abs(tempo[0]-all_atoms[j][0]), 2.0) +
                       Math.pow(Math.abs(tempo[1]-all_atoms[j][1]), 2.0) +
                       Math.pow(Math.abs(tempo[2]-all_atoms[j][2]), 2.0);
            distance = Math.sqrt(distance);
            if (distance < lowest_distance) lowest_distance = distance;

          }
        }
      }
      System.out.printf( "lowest distance = %2.6f \n", lowest_distance );
      if (lowest_distance < 1.0e-4) 
        System.out.printf( "WARNING: The structure has at least 2 atoms that coincide\n");
    
      // Nearest neighbor table
      prn_table.printf("Nearest Neighbor list of neighbors:\n");
      prn_table.printf("atom #    neighbor \n");
      prn_table.printf("------    -------- \n");
      for( i=0; i < natom; i++) {
  
        tempo[0] = all_atoms[i][0];
        tempo[1] = all_atoms[i][1];
        tempo[2] = all_atoms[i][2];
      
        sum_neighbor = 0;
        for( j=0; j < natom; j++) {
          if (j != i) {
            distance = Math.pow(Math.abs(tempo[0]-all_atoms[j][0]), 2.0) +
                       Math.pow(Math.abs(tempo[1]-all_atoms[j][1]), 2.0) +
                       Math.pow(Math.abs(tempo[2]-all_atoms[j][2]), 2.0);
            distance = Math.sqrt(distance);
            if (distance < lowest_distance*(1.4)) {
               table_neighbor[i][sum_neighbor] = j;
               sum_neighbor += 1;
               prn_table.printf("%4d  %4d   (%2.6f)\n", i+1, j+1, distance);
            }
          }
        }
        table_amount[i] = sum_neighbor;
        prn_table.printf("\n");
      }

//    calculation of the number of hydrogen, oxygen, and silicon:
      num_of_hydrogen = 0;
      num_of_oxygen = 0;
      num_of_silicon = 0;
      for( i=0; i < natom; i++) {
        if (table_amount[i] == 1) { num_of_hydrogen += 1; num_of_oxygen  += 1; num_of_silicon -= 1; };
        if (table_amount[i] == 2) {                       num_of_oxygen  += 1; };
        if (table_amount[i] == 3) { num_of_hydrogen += 1; num_of_oxygen  += 1; };
//      if (table_amount[i] == 4) {};
      }
      System.out.printf("num_of_hydrogen=%4d\n num_of_oxygen=%4d\n num_of_silicon=%4d\n natom=%4d\n", 
                         num_of_hydrogen,      num_of_oxygen,      num_of_silicon,      natom);
    
      prn_jmol.printf ("%4d\n\n", natom + num_of_hydrogen + num_of_oxygen + num_of_silicon);
      Data_out.writeInt(natom + num_of_hydrogen + num_of_oxygen + num_of_silicon);

//For the passivation, I assume no symmetry whatsoever except that there are 4 bonds that form a tetrahedra.
  
      for( i=0; i < natom; i++) {
        if (table_amount[i] == 4) {
          prn_jmol.printf ("Si  %.4f  %.4f  %.4f\n",
                     all_atoms[i][0], all_atoms[i][1],all_atoms[i][2]);
          Data_out.writeInt(14);
          Data_out.writeDouble(all_atoms[i][0]);
          Data_out.writeDouble(all_atoms[i][1]);
          Data_out.writeDouble(all_atoms[i][2]);
        }
      }
      for( i=0; i < natom; i++) {
        if (table_amount[i] == 3) {
          prn_jmol.printf ("Si  %.4f  %.4f  %.4f\n",
                     all_atoms[i][0], all_atoms[i][1],all_atoms[i][2]);
          Data_out.writeInt(14);
          Data_out.writeDouble(all_atoms[i][0]);
          Data_out.writeDouble(all_atoms[i][1]);
          Data_out.writeDouble(all_atoms[i][2]);

          j = table_neighbor[i][0];
          k = table_neighbor[i][1];
          m = table_neighbor[i][2];
         
          vectorA[0] = all_atoms[j][0] - all_atoms[k][0];
          vectorA[1] = all_atoms[j][1] - all_atoms[k][1];
          vectorA[2] = all_atoms[j][2] - all_atoms[k][2];
          
          vectorB[0] = all_atoms[m][0] - all_atoms[k][0];
          vectorB[1] = all_atoms[m][1] - all_atoms[k][1];
          vectorB[2] = all_atoms[m][2] - all_atoms[k][2];
          
//        cross product between A and B:
          vectorC[0] = vectorA[1]*vectorB[2] - vectorA[2]*vectorB[1];
          vectorC[1] = vectorA[2]*vectorB[0] - vectorA[0]*vectorB[2];
          vectorC[2] = vectorA[0]*vectorB[1] - vectorA[1]*vectorB[0];
          
          norm = Math.pow(Math.abs(vectorC[0]), 2.0) +
                 Math.pow(Math.abs(vectorC[1]), 2.0) +
                 Math.pow(Math.abs(vectorC[2]), 2.0);
          norm = Math.sqrt(norm);
          vectorC[0] /= norm;
          vectorC[1] /= norm;
          vectorC[2] /= norm;

//        check the sign of the cross product: use the scalar product between the norm and the atom considered 
          vector[0] = all_atoms[i][0] - all_atoms[k][0];
          vector[1] = all_atoms[i][1] - all_atoms[k][1];
          vector[2] = all_atoms[i][2] - all_atoms[k][2];
          
          scalar_product =  vectorC[0]*vector[0] + vectorC[1]*vector[1] + vectorC[2]*vector[2];
          if (scalar_product < 0.0) {
            vectorC[0] *= -1.0;
            vectorC[1] *= -1.0;
            vectorC[2] *= -1.0;
          }

          tempo[0] = all_atoms[i][0] + SiO_bondlength*vectorC[0];
          tempo[1] = all_atoms[i][1] + SiO_bondlength*vectorC[1];
          tempo[2] = all_atoms[i][2] + SiO_bondlength*vectorC[2];

          prn_jmol.printf ("O  %.4f  %.4f  %.4f\n", tempo[0], tempo[1], tempo[2]);
          Data_out.writeInt(8);
          Data_out.writeDouble(tempo[0]);
          Data_out.writeDouble(tempo[1]);
          Data_out.writeDouble(tempo[2]);

          tempo[0] = all_atoms[i][0] + (SiO_bondlength + OH_bondlength)*vectorC[0];
          tempo[1] = all_atoms[i][1] + (SiO_bondlength + OH_bondlength)*vectorC[1];
          tempo[2] = all_atoms[i][2] + (SiO_bondlength + OH_bondlength)*vectorC[2];

          prn_jmol.printf ("H  %.4f  %.4f  %.4f\n", tempo[0], tempo[1], tempo[2]);
          Data_out.writeInt(1);
          Data_out.writeDouble(tempo[0]);
          Data_out.writeDouble(tempo[1]);
          Data_out.writeDouble(tempo[2]);
        }
      }
      for( i=0; i < natom; i++) {
        if (table_amount[i] == 2) {

          prn_jmol.printf ("Si  %.4f  %.4f  %.4f\n",
                     all_atoms[i][0], all_atoms[i][1],all_atoms[i][2]);
          Data_out.writeInt(14);
          Data_out.writeDouble(all_atoms[i][0]);
          Data_out.writeDouble(all_atoms[i][1]);
          Data_out.writeDouble(all_atoms[i][2]);

//        In that case, the 2 hydrogens are replaced by 1 oxygen, much simpler structure in fact

          j = table_neighbor[i][0];
          k = table_neighbor[i][1];

          vectorA[0] = all_atoms[i][0] - all_atoms[j][0];
          vectorA[1] = all_atoms[i][1] - all_atoms[j][1];
          vectorA[2] = all_atoms[i][2] - all_atoms[j][2];
          
          vectorB[0] = all_atoms[i][0] - all_atoms[k][0];
          vectorB[1] = all_atoms[i][1] - all_atoms[k][1];
          vectorB[2] = all_atoms[i][2] - all_atoms[k][2];

          vectorB[0] += vectorA[0];
          vectorB[1] += vectorA[1];
          vectorB[2] += vectorA[2];

          norm = Math.pow(Math.abs(vectorB[0]), 2.0) +
                 Math.pow(Math.abs(vectorB[1]), 2.0) +
                 Math.pow(Math.abs(vectorB[2]), 2.0);
          norm = Math.sqrt(norm);
          vectorB[0] /= norm;
          vectorB[1] /= norm;
          vectorB[2] /= norm;

//        vectorB is parallel to the plane and normalized

          tempo[0] = all_atoms[i][0] + SiO_bondlength*vectorB[0];
          tempo[1] = all_atoms[i][1] + SiO_bondlength*vectorB[1];
          tempo[2] = all_atoms[i][2] + SiO_bondlength*vectorB[2];

          prn_jmol.printf ("O   %.4f  %.4f  %.4f\n", tempo[0], tempo[1], tempo[2]);
          Data_out.writeInt(8);
          Data_out.writeDouble(tempo[0]);
          Data_out.writeDouble(tempo[1]);
          Data_out.writeDouble(tempo[2]);

        }
      }
      for( i=0; i < natom; i++) {
//      In that case, I replace the Silicon atom by a hydrogen and displace the atom at the distance of 1.5 Ang.
        if (table_amount[i] == 1) {
          j = table_neighbor[i][0];

          vector[0] = all_atoms[i][0] - all_atoms[j][0];
          vector[1] = all_atoms[i][1] - all_atoms[j][1];
          vector[2] = all_atoms[i][2] - all_atoms[j][2];
          norm = Math.pow(Math.abs(vector[0]), 2.0) +
                 Math.pow(Math.abs(vector[1]), 2.0) +
                 Math.pow(Math.abs(vector[2]), 2.0);
          norm = Math.sqrt(norm);
          vector[0] /= norm;
          vector[1] /= norm;
          vector[2] /= norm;

          tempo[0] = all_atoms[j][0] + SiO_bondlength*vector[0];
          tempo[1] = all_atoms[j][1] + SiO_bondlength*vector[1];
          tempo[2] = all_atoms[j][2] + SiO_bondlength*vector[2];

          prn_jmol.printf ("O   %.4f  %.4f  %.4f\n", tempo[0], tempo[1], tempo[2]);
          Data_out.writeInt(8);
          Data_out.writeDouble(tempo[0]);
          Data_out.writeDouble(tempo[1]);
          Data_out.writeDouble(tempo[2]);

          tempo[0] = all_atoms[j][0] + (SiO_bondlength + OH_bondlength)*vector[0];
          tempo[1] = all_atoms[j][1] + (SiO_bondlength + OH_bondlength)*vector[1];
          tempo[2] = all_atoms[j][2] + (SiO_bondlength + OH_bondlength)*vector[2];

          prn_jmol.printf ("H   %.4f  %.4f  %.4f\n", tempo[0], tempo[1], tempo[2]);
          Data_out.writeInt(1);
          Data_out.writeDouble(tempo[0]);
          Data_out.writeDouble(tempo[1]);
          Data_out.writeDouble(tempo[2]);
        }
      }

      prn_table.close();
      Data_out.close();

    }
    catch(IOException e){
      System.out.println("Error reading input file");
    }
  }
}