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

    int i, j, natom, ineigh;
    int sum_neighbor;

    FileInputStream fis;
    FileOutputStream Table_out, DataOut_ordered, Jmol_out;

    PrintStream prn_table, prn_jmol; 

    int [] atom_number;
    double [][] all_atoms;
    int [][] table_neighbor;
    int [] table_amount;
    double lowest_distanceAtom = 1.0e+20;
    double lowest_distanceHy = 1.0e+20;
    double distanceAtom;
    double distanceHy;
    double [] tempo;
    double au = 1.0; // 0.529177;


    tempo = new double [3];

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

      Table_out = new FileOutputStream("NNtable_WITH_hydrogen.list");
      Jmol_out = new FileOutputStream("ordered_WITH_hydrogen.for_parsec");
      DataOut_ordered = new FileOutputStream("ordered_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_number   = new int [natom];
      table_amount   = new int [natom];
      table_neighbor = new int [natom][6];

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

      // find lowest distances
      for( i=0; i < natom; i++) {
          tempo[0] = all_atoms[i][0];
          tempo[1] = all_atoms[i][1];
          tempo[2] = all_atoms[i][2];
        if (atom_number[i] != 1) {

          for( j=0; j < natom; j++) {
            if (j != i && atom_number[j] != 1) {
              distanceAtom = 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);
              distanceAtom = Math.sqrt(distanceAtom);
              if (distanceAtom < lowest_distanceAtom) lowest_distanceAtom = distanceAtom;
            }       
          }
        }
        if (atom_number[i] == 1) {

          for( j=0; j < natom; j++) {
            if (j != i ) {
              distanceHy = 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);
              distanceHy = Math.sqrt(distanceHy);
              if (distanceHy < lowest_distanceHy) lowest_distanceHy = distanceHy;
            }       
          }
        }
      }
      System.out.printf( "lowest distance Si= %2.6f \n", lowest_distanceAtom );
      System.out.printf( "lowest distance Hy= %2.6f \n", lowest_distanceHy );
      if (lowest_distanceAtom < 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];
      
        if (atom_number[i] != 1) {
          sum_neighbor = 0;
          for( j=0; j < natom; j++) {
            if (j != i) {
              distanceAtom = 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);
              distanceAtom = Math.sqrt(distanceAtom);
              if (distanceAtom < lowest_distanceAtom*(1.05)) {
                 table_neighbor[i][sum_neighbor] = j;
                 sum_neighbor += 1;
                 prn_table.printf("%4d  %4d   (%2.6f)\n", i+1, j+1, distanceAtom);
              }
            }
          }
          table_amount[i] = sum_neighbor;
        } else {
          sum_neighbor = 0;
          for( j=0; j < natom; j++) {
            if (j != i) {
              distanceHy = 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);
              distanceHy = Math.sqrt(distanceHy);
              if (distanceHy < lowest_distanceHy*(1.05)) {
                 table_neighbor[i][sum_neighbor] = j;
                 sum_neighbor += 1;
                 prn_table.printf("%4d  %4d   (%2.6f)\n", i+1, j+1, distanceHy);
              }
            }
          }
          table_amount[i] = sum_neighbor;
        }
        prn_table.printf("\n");
      }

      prn_jmol.printf ("%4d\n\n", natom);
      prn_jmol.printf("Lead\n");
      for( i=0; i < natom; i++) {
        if (atom_number[i] == 82) {
          prn_table.printf("%4d: %4d   %4d   %4d   %4d\n", i+1, table_neighbor[i][0]+1, 
                                                                table_neighbor[i][1]+1, 
                                                                table_neighbor[i][2]+1, 
                                                                table_neighbor[i][3]+1);
          prn_jmol.printf ("Pb  %3.6f  %3.6f  %3.6f\n",
                   all_atoms[i][0]/au, all_atoms[i][1]/au,all_atoms[i][2]/au);
        }
      }
      prn_jmol.printf("Selenium\n");
      for( i=0; i < natom; i++) {
        if (atom_number[i] == 34) {
          prn_table.printf("%4d: %4d   %4d   %4d   %4d\n", i+1, table_neighbor[i][0]+1, 
                                                                table_neighbor[i][1]+1, 
                                                                table_neighbor[i][2]+1, 
                                                                table_neighbor[i][3]+1);
          prn_jmol.printf ("Se  %3.6f  %3.6f  %3.6f\n",
                   all_atoms[i][0]/au, all_atoms[i][1]/au,all_atoms[i][2]/au);
        }
      }
      prn_jmol.printf("Hydrogen bonded to Cadmium\n");
      for( i=0; i < natom; i++) {
        if (atom_number[i] == 1) {
          prn_table.printf("%4d: %4d\n", i+1, table_neighbor[i][0]+1);
          ineigh = table_neighbor[i][0];
          if (atom_number[ineigh] == 82) {
            prn_jmol.printf ("H  %3.6f  %3.6f  %3.6f\n",
                             all_atoms[i][0]/au, all_atoms[i][1]/au,all_atoms[i][2]/au);
          }
        }
      }
      prn_jmol.printf("Hydrogen bonded to Selenium\n");
      for( i=0; i < natom; i++) {
        if (atom_number[i] == 1) {
          prn_table.printf("%4d: %4d\n", i+1, table_neighbor[i][0]+1);
          ineigh = table_neighbor[i][0];
          if (atom_number[ineigh] == 34) {
            prn_jmol.printf ("He %3.6f  %3.6f  %3.6f\n",
                             all_atoms[i][0]/au, all_atoms[i][1]/au,all_atoms[i][2]/au);
          }
        }
      }

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

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