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

    char aa;
    int i, j, k, m, natom, neighbor, center_atom;
    int atom1, atom2, parallel, atom_dangling;
    int dangling, other_neighbor;
    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_atoms;
    double lowest_distance = 1.0e+20;
    double distance, norm, scalar_product;
    double [] tempo, vector, vectorA, vectorB, vectorC;

    double SiH_bondlength = 1.1;

    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_hydrogen.list");
      Jmol_out = new FileOutputStream("hydrogen_passivated.xyz");
      DataOut_ordered = new FileOutputStream("WITH_hydrogen.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];
      atom_name = new String [natom];
      table_amount   = new int [natom];
      table_neighbor = new int [natom][6];

      for( i=0; i < natom; i++) {
          atom_name[i] = 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.05)) {
               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 hydrogens and atoms (cations + anions):
      num_of_hydrogen = 0;
      num_of_atoms = 0;
      for( i=0; i < natom; i++) {
        if (table_amount[i] == 1) { num_of_hydrogen += 1; num_of_atoms -= 1; }
        if (table_amount[i] == 2) { num_of_hydrogen += 2; num_of_atoms -= 1; }
        if (table_amount[i] == 3)   num_of_hydrogen += 3;
        if (table_amount[i] == 4)   num_of_hydrogen += 2;
        if (table_amount[i] == 5)   num_of_hydrogen += 1;
      }
      System.out.printf("num_of_hydrogen=%4d, num_of_atoms=%4d, natom=%4d\n", num_of_hydrogen, num_of_atoms, natom);
    
      prn_jmol.printf ("%4d\n\n", natom + num_of_hydrogen + num_of_atoms);
      Data_out.writeInt(natom + num_of_hydrogen + num_of_atoms);

//For the passivation, I assume no symmetry whatsoever except that there are 6 possible bonds.
  
      for( center_atom=0; center_atom < natom; center_atom++) {
        if (table_amount[center_atom] == 6) {
          // no passivation necessary in this case
          aa = atom_name[center_atom].charAt(0);
          if (aa == 'P') {
            prn_jmol.printf ("Pb  %.4f  %.4f  %.4f\n",
                       all_atoms[center_atom][0], all_atoms[center_atom][1],all_atoms[center_atom][2]);
            Data_out.writeInt(82);
          };
          if (aa == 'S') {
            prn_jmol.printf ("Se  %.4f  %.4f  %.4f\n",
                       all_atoms[center_atom][0], all_atoms[center_atom][1],all_atoms[center_atom][2]);
            Data_out.writeInt(34);
          };
          Data_out.writeDouble(all_atoms[center_atom][0]);
          Data_out.writeDouble(all_atoms[center_atom][1]);
          Data_out.writeDouble(all_atoms[center_atom][2]);
        }
      }

      System.out.printf( "1 dangling bond\n" );
      for( center_atom=0; center_atom < natom; center_atom++) {
        if (table_amount[center_atom] == 5) {
          aa = atom_name[center_atom].charAt(0);
          if (aa == 'P') {
            prn_jmol.printf ("Pb  %.4f  %.4f  %.4f\n",
                       all_atoms[center_atom][0], all_atoms[center_atom][1],all_atoms[center_atom][2]);
            Data_out.writeInt(82);
          };
          if (aa == 'S') {
            prn_jmol.printf ("Se  %.4f  %.4f  %.4f\n",
                       all_atoms[center_atom][0], all_atoms[center_atom][1],all_atoms[center_atom][2]);
            Data_out.writeInt(34);
          };
          Data_out.writeDouble(all_atoms[center_atom][0]);
          Data_out.writeDouble(all_atoms[center_atom][1]);
          Data_out.writeDouble(all_atoms[center_atom][2]);

//        The sole dangling bond is opposite to one atom orthogonal to all the rest of the atoms
//        We use that as a condition for finding the dangling bond:
//        Do a loop on all atoms, one after the other;
//        if all atom are orthogonal, then that atom is opposite to the dangling bond, simply.
         
          //Choose a first atom
          atom_dangling = table_neighbor[center_atom][0];
          for (neighbor = 0; neighbor < 5; neighbor++) {
             parallel = 0;
             atom1 = table_neighbor[center_atom][neighbor];
             vectorA[0] = all_atoms[atom1][0] - all_atoms[center_atom][0];
             vectorA[1] = all_atoms[atom1][1] - all_atoms[center_atom][1];
             vectorA[2] = all_atoms[atom1][2] - all_atoms[center_atom][2];
          
             //verify scalar product with first atom
             for (other_neighbor = 0; other_neighbor < 5; other_neighbor++) {
                if (other_neighbor != neighbor) {
                   atom2 = table_neighbor[center_atom][other_neighbor];
                   vectorB[0] = all_atoms[atom2][0] - all_atoms[center_atom][0];
                   vectorB[1] = all_atoms[atom2][1] - all_atoms[center_atom][1];
                   vectorB[2] = all_atoms[atom2][2] - all_atoms[center_atom][2];
                
                   norm =  vectorA[0]*vectorB[0] + vectorA[1]*vectorB[1] + vectorA[2]*vectorB[2];
                   if (Math.abs(norm) > 0.5) {
                      parallel = 1;
                      break; // found that it's not it; stop this loop
                   }
                }
             }
             if (parallel == 0) {
                dangling = neighbor;
                atom_dangling = table_neighbor[center_atom][dangling];
                break; // found it; stop this loop
             }
          }
 
          System.out.printf( "atom_dangling = %4d center_atom = %4d\n", atom_dangling+1, center_atom+1 );
          // put an hydrogen opposite to the dangling bond atom found
          vectorA[0] = -all_atoms[atom_dangling][0] + all_atoms[center_atom][0];
          vectorA[1] = -all_atoms[atom_dangling][1] + all_atoms[center_atom][1];
          vectorA[2] = -all_atoms[atom_dangling][2] + all_atoms[center_atom][2];
          norm = Math.pow(Math.abs(vectorA[0]), 2.0) +
                 Math.pow(Math.abs(vectorA[1]), 2.0) +
                 Math.pow(Math.abs(vectorA[2]), 2.0);
          norm = Math.sqrt(norm);
          vectorA[0] /= norm;
          vectorA[1] /= norm;
          vectorA[2] /= norm;

          tempo[0] = all_atoms[center_atom][0] + SiH_bondlength*vectorA[0];
          tempo[1] = all_atoms[center_atom][1] + SiH_bondlength*vectorA[1];
          tempo[2] = all_atoms[center_atom][2] + SiH_bondlength*vectorA[2];


          if (aa == 'P') {
            prn_jmol.printf ("H  %.4f  %.4f  %.4f\n", tempo[0], tempo[1], tempo[2]);
            Data_out.writeInt(1);
          };
          if (aa == 'S') {
            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]);
        }
      }

      System.out.printf( "2 dangling bond\n" );
      for( center_atom=0; center_atom < natom; center_atom++) {
        if (table_amount[center_atom] == 4) {
          aa = atom_name[center_atom].charAt(0);
          if (aa == 'P') {
            prn_jmol.printf ("Pb  %.4f  %.4f  %.4f\n",
                       all_atoms[center_atom][0], all_atoms[center_atom][1],all_atoms[center_atom][2]);
            Data_out.writeInt(82);
          };
          if (aa == 'S') {
            prn_jmol.printf ("Se  %.4f  %.4f  %.4f\n",
                       all_atoms[center_atom][0], all_atoms[center_atom][1],all_atoms[center_atom][2]);
            Data_out.writeInt(34);
          };
          Data_out.writeDouble(all_atoms[center_atom][0]);
          Data_out.writeDouble(all_atoms[center_atom][1]);
          Data_out.writeDouble(all_atoms[center_atom][2]);

//        Whenever I find an atom that has no opposite I put an atom;
//        Equivalently, for all atoms that has an opposite, skip.
//        That means when scalar product is < 0 (-1*norm)

          //Initialize a first atom that is just to have java to compile; compiler more careful than me!
          atom_dangling = table_neighbor[center_atom][0];
          for (neighbor = 0; neighbor < 4; neighbor++) {
             parallel = 0;
             atom1 = table_neighbor[center_atom][neighbor];
             vectorA[0] = all_atoms[atom1][0] - all_atoms[center_atom][0];
             vectorA[1] = all_atoms[atom1][1] - all_atoms[center_atom][1];
             vectorA[2] = all_atoms[atom1][2] - all_atoms[center_atom][2];
          
             //verify scalar product with that atom in the loop
             for (other_neighbor = 0; other_neighbor < 4; other_neighbor++) {
                if (other_neighbor != neighbor) {
                   atom2 = table_neighbor[center_atom][other_neighbor];
                   vectorB[0] = all_atoms[atom2][0] - all_atoms[center_atom][0];
                   vectorB[1] = all_atoms[atom2][1] - all_atoms[center_atom][1];
                   vectorB[2] = all_atoms[atom2][2] - all_atoms[center_atom][2];
                
                   norm =  vectorA[0]*vectorB[0] + vectorA[1]*vectorB[1] + vectorA[2]*vectorB[2];
                   if (Math.abs(norm) > 0.5) {
                      parallel = 1;
                      break; // found that it's not it; stop this loop
                   }
                }
             }
             if (parallel == 0) {
                dangling = neighbor;
                atom_dangling = table_neighbor[center_atom][dangling];
                // the difference with above is that here I don't break
                // that means that we loop on all the atoms
                System.out.printf( "atom_dangling = %4d center_atom = %4d\n", atom_dangling+1, center_atom+1 );
                // put an hydrogen opposite to the dangling bond atom found
                vectorA[0] = -all_atoms[atom_dangling][0] + all_atoms[center_atom][0];
                vectorA[1] = -all_atoms[atom_dangling][1] + all_atoms[center_atom][1];
                vectorA[2] = -all_atoms[atom_dangling][2] + all_atoms[center_atom][2];
                norm = Math.pow(Math.abs(vectorA[0]), 2.0) +
                       Math.pow(Math.abs(vectorA[1]), 2.0) +
                       Math.pow(Math.abs(vectorA[2]), 2.0);
                norm = Math.sqrt(norm);
                vectorA[0] /= norm;
                vectorA[1] /= norm;
                vectorA[2] /= norm;
      
                tempo[0] = all_atoms[center_atom][0] + SiH_bondlength*vectorA[0];
                tempo[1] = all_atoms[center_atom][1] + SiH_bondlength*vectorA[1];
                tempo[2] = all_atoms[center_atom][2] + SiH_bondlength*vectorA[2];
      
      
                if (aa == 'P') {
                  prn_jmol.printf ("H  %.4f  %.4f  %.4f\n", tempo[0], tempo[1], tempo[2]);
                  Data_out.writeInt(1);
                };
                if (aa == 'S') {
                  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]);
             }
          }
        }
      }

      System.out.printf( "3 dangling bond\n" );
      for( center_atom=0; center_atom < natom; center_atom++) {
        if (table_amount[center_atom] == 3) {
          aa = atom_name[center_atom].charAt(0);
          if (aa == 'P') {
            prn_jmol.printf ("Pb  %.4f  %.4f  %.4f\n",
                       all_atoms[center_atom][0], all_atoms[center_atom][1],all_atoms[center_atom][2]);
            Data_out.writeInt(82);
          };
          if (aa == 'S') {
            prn_jmol.printf ("Se  %.4f  %.4f  %.4f\n",
                       all_atoms[center_atom][0], all_atoms[center_atom][1],all_atoms[center_atom][2]);
            Data_out.writeInt(34);
          };
          Data_out.writeDouble(all_atoms[center_atom][0]);
          Data_out.writeDouble(all_atoms[center_atom][1]);
          Data_out.writeDouble(all_atoms[center_atom][2]);

//        Whenever I find an atom that has no opposite I put an atom;
//        Equivalently, for all atoms that has an opposite, skip.
//        That means when scalar product is < 0 (-1*norm)

          //Initialize a first atom that is just to have java to compile; compiler more careful than me!
          atom_dangling = table_neighbor[center_atom][0];
          for (neighbor = 0; neighbor < 3; neighbor++) {
             parallel = 0;
             atom1 = table_neighbor[center_atom][neighbor];
             vectorA[0] = all_atoms[atom1][0] - all_atoms[center_atom][0];
             vectorA[1] = all_atoms[atom1][1] - all_atoms[center_atom][1];
             vectorA[2] = all_atoms[atom1][2] - all_atoms[center_atom][2];
          
             //verify scalar product with that atom in the loop
             for (other_neighbor = 0; other_neighbor < 3; other_neighbor++) {
                if (other_neighbor != neighbor) {
                   atom2 = table_neighbor[center_atom][other_neighbor];
                   vectorB[0] = all_atoms[atom2][0] - all_atoms[center_atom][0];
                   vectorB[1] = all_atoms[atom2][1] - all_atoms[center_atom][1];
                   vectorB[2] = all_atoms[atom2][2] - all_atoms[center_atom][2];
                
                   norm =  vectorA[0]*vectorB[0] + vectorA[1]*vectorB[1] + vectorA[2]*vectorB[2];
                   if (Math.abs(norm) > 0.5) {
                      parallel = 1;
                      break; // found that it's not it; stop this loop
                   }
                }
             }
             if (parallel == 0) {
                dangling = neighbor;
                atom_dangling = table_neighbor[center_atom][dangling];
                // the difference with above is that here I don't break
                // that means that we loop on all the atoms
                System.out.printf( "atom_dangling = %4d center_atom = %4d\n", atom_dangling+1, center_atom+1 );
                // put an hydrogen opposite to the dangling bond atom found
                vectorA[0] = -all_atoms[atom_dangling][0] + all_atoms[center_atom][0];
                vectorA[1] = -all_atoms[atom_dangling][1] + all_atoms[center_atom][1];
                vectorA[2] = -all_atoms[atom_dangling][2] + all_atoms[center_atom][2];
                norm = Math.pow(Math.abs(vectorA[0]), 2.0) +
                       Math.pow(Math.abs(vectorA[1]), 2.0) +
                       Math.pow(Math.abs(vectorA[2]), 2.0);
                norm = Math.sqrt(norm);
                vectorA[0] /= norm;
                vectorA[1] /= norm;
                vectorA[2] /= norm;
      
                tempo[0] = all_atoms[center_atom][0] + SiH_bondlength*vectorA[0];
                tempo[1] = all_atoms[center_atom][1] + SiH_bondlength*vectorA[1];
                tempo[2] = all_atoms[center_atom][2] + SiH_bondlength*vectorA[2];
      
      
                if (aa == 'P') {
                  prn_jmol.printf ("H  %.4f  %.4f  %.4f\n", tempo[0], tempo[1], tempo[2]);
                  Data_out.writeInt(1);
                };
                if (aa == 'S') {
                  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]);
             }
          }

        }
      }

      System.out.printf( "4 dangling bond\n" );
      for( center_atom=0; center_atom < natom; center_atom++) {
        if (table_amount[center_atom] == 2) {
          aa = atom_name[center_atom].charAt(0);
//        if (aa == 'P') {
//          prn_jmol.printf ("Pb  %.4f  %.4f  %.4f\n",
//                     all_atoms[center_atom][0], all_atoms[center_atom][1],all_atoms[center_atom][2]);
//          Data_out.writeInt(82);
//        };
//        if (aa == 'S') {
//          prn_jmol.printf ("Se  %.4f  %.4f  %.4f\n",
//                     all_atoms[center_atom][0], all_atoms[center_atom][1],all_atoms[center_atom][2]);
//          Data_out.writeInt(34);
//        };
//        Data_out.writeDouble(all_atoms[center_atom][0]);
//        Data_out.writeDouble(all_atoms[center_atom][1]);
//        Data_out.writeDouble(all_atoms[center_atom][2]);

//        Whenever I find an atom that has no opposite I put an atom;
//        Equivalently, for all atoms that has an opposite, skip.
//        That means when scalar product is < 0 (-1*norm)

          //Initialize a first atom that is just to have java to compile; compiler more careful than me!
          atom_dangling = table_neighbor[center_atom][0];
          for (neighbor = 0; neighbor < 2; neighbor++) {
             parallel = 0;
             atom1 = table_neighbor[center_atom][neighbor];
             vectorA[0] = all_atoms[atom1][0] - all_atoms[center_atom][0];
             vectorA[1] = all_atoms[atom1][1] - all_atoms[center_atom][1];
             vectorA[2] = all_atoms[atom1][2] - all_atoms[center_atom][2];
          
             //verify scalar product with that atom in the loop
             for (other_neighbor = 0; other_neighbor < 2; other_neighbor++) {
                if (other_neighbor != neighbor) {
                   atom2 = table_neighbor[center_atom][other_neighbor];
                   vectorB[0] = all_atoms[atom2][0] - all_atoms[center_atom][0];
                   vectorB[1] = all_atoms[atom2][1] - all_atoms[center_atom][1];
                   vectorB[2] = all_atoms[atom2][2] - all_atoms[center_atom][2];
                
                   norm =  vectorA[0]*vectorB[0] + vectorA[1]*vectorB[1] + vectorA[2]*vectorB[2];
                   if (Math.abs(norm) > 0.5) {
                      parallel = 1;
                      break; // found that it's not it; stop this loop
                   }
                }
             }
             if (parallel == 0) {
                dangling = neighbor;
                atom_dangling = table_neighbor[center_atom][dangling];
                // the difference with above is that here I don't break
                // that means that we loop on all the atoms
                System.out.printf( "atom_dangling = %4d center_atom = %4d\n", atom_dangling+1, center_atom+1 );
                // put an hydrogen opposite to the dangling bond atom found
                vectorA[0] = -all_atoms[atom_dangling][0] + all_atoms[center_atom][0];
                vectorA[1] = -all_atoms[atom_dangling][1] + all_atoms[center_atom][1];
                vectorA[2] = -all_atoms[atom_dangling][2] + all_atoms[center_atom][2];
                norm = Math.pow(Math.abs(vectorA[0]), 2.0) +
                       Math.pow(Math.abs(vectorA[1]), 2.0) +
                       Math.pow(Math.abs(vectorA[2]), 2.0);
                norm = Math.sqrt(norm);
                vectorA[0] /= norm;
                vectorA[1] /= norm;
                vectorA[2] /= norm;
      
                tempo[0] = all_atoms[atom_dangling][0] + SiH_bondlength*vectorA[0];
                tempo[1] = all_atoms[atom_dangling][1] + SiH_bondlength*vectorA[1];
                tempo[2] = all_atoms[atom_dangling][2] + SiH_bondlength*vectorA[2];
      
      
                if (aa == 'P') {
                  prn_jmol.printf ("H  %.4f  %.4f  %.4f\n", tempo[0], tempo[1], tempo[2]);
                  Data_out.writeInt(1);
                };
                if (aa == 'S') {
                  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]);
             }
          }

        }
      }

      System.out.printf( "5 dangling bond\n" );
      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];
          aa = atom_name[i].charAt(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] + SiH_bondlength*vector[0];
          tempo[1] = all_atoms[j][1] + SiH_bondlength*vector[1];
          tempo[2] = all_atoms[j][2] + SiH_bondlength*vector[2];

          if (aa == 'P') {
            prn_jmol.printf ("H  %.4f  %.4f  %.4f\n", tempo[0], tempo[1], tempo[2]);
            Data_out.writeInt(1);
          };
          if (aa == 'S') {
            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");
    }
  }
}