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

    int i, j, natom;
    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_distanceSi = 1.0e+20;
    double lowest_distanceOx = 1.0e+20;
    double lowest_distanceHy = 1.0e+20;
    double distanceSi;
    double distanceOx;
    double distanceHy;
    double [] tempo;


    tempo = new double [3];

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

      Table_out = new FileOutputStream("NNtable_WITH_oxygen.list");
      Jmol_out = new FileOutputStream("ordered_WITH_oxygen.for_parsec");
      DataOut_ordered = new FileOutputStream("ordered_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];
      atom_number   = new int [natom];
      table_amount   = new int [natom];
      table_neighbor = new int [natom][4];

      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();
//        System.out.printf("%4d    %2.6f   %2.6f   %2.6f\n", atom_number[i], all_atoms[i][0], all_atoms[i][1], all_atoms[i][2]);
      }

      // 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] == 14) {

            for( j=0; j < natom; j++) {
              if (j != i && atom_number[j] == 14) {
                distanceSi = 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);
                distanceSi = Math.sqrt(distanceSi);
                if (distanceSi < lowest_distanceSi) lowest_distanceSi = distanceSi;
              }       
              if (j != i && atom_number[j] == 8) {
                distanceOx = 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);
                distanceOx = Math.sqrt(distanceOx);
                if (distanceOx < lowest_distanceOx) lowest_distanceOx = distanceOx;
              }       
            }
          }
          if (atom_number[i] == 8) {
            for( j=0; j < natom; j++) {
              if (j != i && atom_number[j] == 1) {
                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-Si= %2.6f \n", lowest_distanceSi );
      System.out.printf( "lowest distance O - O= %2.6f \n", lowest_distanceOx );
      System.out.printf( "lowest distance O - H= %2.6f \n", lowest_distanceHy );
      if (lowest_distanceSi < 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] == 14) {
          sum_neighbor = 0;
          for( j=0; j < natom; j++) {
            if (j != i) {
              distanceSi = 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);
              distanceSi = Math.sqrt(distanceSi);
              if (distanceSi < lowest_distanceSi*(1.1)) {
                 table_neighbor[i][sum_neighbor] = j;
                 sum_neighbor += 1;
                 prn_table.printf("%4d  %4d   (%2.6f)\n", i+1, j+1, distanceSi);
              }
            }
          }
          table_amount[i] = sum_neighbor;
        }
        if (atom_number[i] == 8) {
          sum_neighbor = 0;
          for( j=0; j < natom; j++) {
            if (j != i) {
              distanceOx = 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);
              distanceOx = Math.sqrt(distanceOx);
              if (distanceOx < lowest_distanceOx*(1.1)) {
                 table_neighbor[i][sum_neighbor] = j;
                 sum_neighbor += 1;
                 prn_table.printf("%4d  %4d   (%2.6f)\n", i+1, j+1, distanceOx);
              }
            }
          }
          table_amount[i] = sum_neighbor;
        }
        if (atom_number[i] == 1) {
          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.1)) {
                 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 ("Silicon\n");
      for( i=0; i < natom; i++) {
        if (atom_number[i] == 14) {
          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 ("%3.6f  %3.6f  %3.6f\n",
                   all_atoms[i][0], all_atoms[i][1],all_atoms[i][2]);
        }
      }
      prn_jmol.printf ("\nOxygen\n");
      for( i=0; i < natom; i++) {
        if (atom_number[i] == 8) {
          prn_table.printf("%4d: %4d   %4d\n", i+1, table_neighbor[i][0]+1,
                                                                table_neighbor[i][1]+1);
          prn_jmol.printf ("%3.6f  %3.6f  %3.6f\n",
                   all_atoms[i][0], all_atoms[i][1],all_atoms[i][2]);
        }
      }
      prn_jmol.printf ("\nHydrogen\n");
      for( i=0; i < natom; i++) {
        if (atom_number[i] == 1) {
          prn_table.printf("%4d: %4d   %4d\n", i+1, table_neighbor[i][0]+1,
                                                                table_neighbor[i][1]+1);
          prn_jmol.printf ("%3.6f  %3.6f  %3.6f\n",
                   all_atoms[i][0], all_atoms[i][1],all_atoms[i][2]);
        }
      }



      Data_out.close();

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