import java.io.*;
import java.lang.Math;

public class GrapheneOnQuartz2jmol {
  public static void main( String[] args ) {
    int i,j,k, icoord, iatom;

    int nXcell = 2, nYcell = 2, nZcell = 1;
    double aa = 4.9134, cc = 5.4052;

    double x,y,z, x_cart, y_cart, z_cart;
    double [][] base, silicon, oxygene, Si_cart, Ox_cart; 
    double [][] carbone, C_cart;
    double []   vector, tempo;
                FileOutputStream out; // declare a file output object
                FileOutputStream parsec_out; // declare a file output object
                PrintStream p; // declare a print stream object

    System.out.println("This program generates the structure of alpha-Quartz");
    System.out.println("Taking the parameters from the Navy");
    System.out.println("Written by carrier@cs.umn.edu");
    System.out.println("Distance C-O= 1.43 Ang");
    System.out.println("Group #154; P6(_3)/mmc");

    int num_atoms_C  = 8;
    int num_atoms_Si = 3;
    int num_atoms_Ox = 6;

//  There are several choices of shifts, these 2 are the most symmetric ones
//  It leads to 2 types of rings
//  double shiftX =  0.16666666;
//  double shiftY = -0.16666667;

    double shiftX =  0.0;
    double shiftY =  0.0;

    carbone = new double [num_atoms_C ][3];
    silicon = new double [num_atoms_Si][3];
    oxygene = new double [num_atoms_Ox][3];

    C_cart  = new double [num_atoms_C ][3];
    Si_cart = new double [num_atoms_Si][3];
    Ox_cart = new double [num_atoms_Ox][3];

    base    = new double [3][3];
    vector  = new double [3];
    tempo   = new double [3];

// The basis
   base[0][0] = 0.5; 
   base[0][1] = -Math.sqrt(3.0)/2.0; 
   base[0][2] = 0.000;

   base[1][0] = 0.5; 
   base[1][1] =  Math.sqrt(3.0)/2.0; 
   base[1][2] = 0.000;

   base[2][0] = 0.0; 
   base[2][1] = 0.0; 
   base[2][2] = cc/aa;
    

// print the basis
    System.out.println("The lattice basis");
    System.out.printf( "%2.12f %2.12f %2.12f\n",base[0][0], base[0][1], base[0][2] );
    System.out.printf( "%2.12f %2.12f %2.12f\n",base[1][0], base[1][1], base[1][2] );
    System.out.printf( "%2.12f %2.12f %2.12f\n",base[2][0], base[2][1], base[2][2] );

// prints to file out
    try
    {
      // Create a new file output stream
      out = new FileOutputStream("GrapheneOnQuartz.xyz");

      // Connect print stream to the output stream
      p = new PrintStream( out );

      p.println(num_atoms_C + num_atoms_Si + num_atoms_Ox);
      p.println();
 
//    makes a supercell of twice the lattice constant of quartz
//    shift = the amount of shift in x and y
//    This number will vary from models to models
      carbone[0][0] = (0.            + shiftX); carbone[0][1] = (0.            + shiftY); carbone[0][2] = (1.05);
      carbone[1][0] = (0.    + 1./2. + shiftX); carbone[1][1] = (0.            + shiftY); carbone[1][2] = (1.05);
      carbone[2][0] = (0.            + shiftX); carbone[2][1] = (0.    + 1./2. + shiftY); carbone[2][2] = (1.05);
      carbone[3][0] = (0.    + 1./2. + shiftX); carbone[3][1] = (0.    + 1./2. + shiftY); carbone[3][2] = (1.05);

      carbone[4][0] = (1./6.         + shiftX); carbone[4][1] = (2./6.         + shiftY); carbone[4][2] = (1.05);
      carbone[5][0] = (1./6. + 1./2. + shiftX); carbone[5][1] = (2./6.         + shiftY); carbone[5][2] = (1.05);
      carbone[6][0] = (1./6.         + shiftX); carbone[6][1] = (2./6. + 1./2. + shiftY); carbone[6][2] = (1.05);
      carbone[7][0] = (1./6. + 1./2. + shiftX); carbone[7][1] = (2./6. + 1./2. + shiftY); carbone[7][2] = (1.05);

      for( i = 0; i < num_atoms_C; i++ ) {
         C_cart[i][0] = aa*(carbone[i][0]*base[0][0] + carbone[i][1]*base[1][0] + carbone[i][2]*base[2][0]);
         C_cart[i][1] = aa*(carbone[i][0]*base[0][1] + carbone[i][1]*base[1][1] + carbone[i][2]*base[2][1]);
         C_cart[i][2] = aa*(carbone[i][0]*base[0][2] + carbone[i][1]*base[1][2] + carbone[i][2]*base[2][2]);
      }

      for( i = 0; i< num_atoms_C; i++ ) {
         p.printf ("C  %2.8f  %2.8f  %2.8f\n",C_cart[i][0], C_cart[i][1], C_cart[i][2]);
      }

      System.out.println( "6c ; O" );
//    Data from  Mineralogy web page
//    x = 0.4135;
//    y = 0.2669;
//    z = 0.1191; // That one makes no sense as it is. According to the military, it is 0.7854 which is +2/3 
                  // the value given by Levien-Prewitt

//    Data from the Navy 
      x = 0.4141;
      y = 0.2681;
      z = 0.1187 + 0.6666;  // That one does fit with the definition of Levien-Prewitt if 2/3 is added

//    Int. Tables for Cryst. #154, 6c
      oxygene[ 0][0] = (    x        ); oxygene[ 0][1] = (      y      ); oxygene[ 0][2] = (    z            ); // 09
      oxygene[ 1][0] = (     -y      ); oxygene[ 1][1] = (    x-y      ); oxygene[ 1][2] = (    z+2./3. - 1.0); // 10
      oxygene[ 2][0] = (   -x+y + 1.0); oxygene[ 2][1] = (   -x   + 1.0); oxygene[ 2][2] = (    z+1./3. - 1.0); // 11
      oxygene[ 3][0] = (      y      ); oxygene[ 3][1] = (    x        ); oxygene[ 3][2] = (   -z       + 1.0); // 12
      oxygene[ 4][0] = (   -x   + 1.0); oxygene[ 4][1] = (   -x+y + 1.0); oxygene[ 4][2] = (   -z+2./3. + 0.0); // 13
      oxygene[ 5][0] = (    x-y      ); oxygene[ 5][1] = (     -y      ); oxygene[ 5][2] = (   -z+1./3. + 1.0); // 14

//    Check for the elements outside the unit cell
//    for( i = 0; i < num_atoms_Ox; i++ ) {
//       for( icoord=0; icoord <=2; icoord++) {
//          if (oxygene[i][icoord] >=  1.0) oxygene[i][icoord] -= 1.0;
//          if (oxygene[i][icoord] <= -1.0) oxygene[i][icoord] += 1.0;
//       }
//    }
//    Transform in cartesian coordinates
      for( i = 0; i < num_atoms_Ox; i++ ) {
         Ox_cart[i][0] = aa*(oxygene[i][0]*base[0][0] + oxygene[i][1]*base[1][0] + oxygene[i][2]*base[2][0]);
         Ox_cart[i][1] = aa*(oxygene[i][0]*base[0][1] + oxygene[i][1]*base[1][1] + oxygene[i][2]*base[2][1]);
         Ox_cart[i][2] = aa*(oxygene[i][0]*base[0][2] + oxygene[i][1]*base[1][2] + oxygene[i][2]*base[2][2]);
      }

      for( i = 0; i < num_atoms_Ox; i++ ) {
         p.printf ("O   %2.8f  %2.8f  %2.8f\n", Ox_cart[i][0], Ox_cart[i][1], Ox_cart[i][2]);
      }

      System.out.println( "3 a ; Si" );
      x = 0.4699 ;

//    Int. Tables for Cryst. #154, 3a:
      silicon[0][0] = (    x        ); silicon[0][1] = (    0        ); silicon[0][2] = (    2./3.    );
      silicon[1][0] = (    0        ); silicon[1][1] = (    x        ); silicon[1][2] = (    1./3.    );
      silicon[2][0] = (   -x + 1.0  ); silicon[2][1] = (   -x + 1.0  ); silicon[2][2] = (    0.       );

//    Check for the elements outside the unit cell
//    for( i = 0; i < num_atoms_Si; i++ ) {
//       for( icoord=0; icoord <=2; icoord++) {
//          if (silicon[i][icoord] >=  1.0) silicon[i][icoord] -= 1.0;
//          if (silicon[i][icoord] <= -1.0) silicon[i][icoord] += 1.0;
//       }
//    }
//    Transform in cartesian coordinates
      for( i = 0; i < num_atoms_Si; i++ ) {
         Si_cart[i][0] = aa*(silicon[i][0]*base[0][0] + silicon[i][1]*base[1][0] + silicon[i][2]*base[2][0]);
         Si_cart[i][1] = aa*(silicon[i][0]*base[0][1] + silicon[i][1]*base[1][1] + silicon[i][2]*base[2][1]);
         Si_cart[i][2] = aa*(silicon[i][0]*base[0][2] + silicon[i][1]*base[1][2] + silicon[i][2]*base[2][2]);
      }

      for( i = 0; i< num_atoms_Si; i++ ) {
         p.printf ("Si   %2.8f  %2.8f  %2.8f\n", Si_cart[i][0], Si_cart[i][1], Si_cart[i][2]);
      }

      p.close();
    }
    catch (Exception e)
    {
      System.err.println ("Error writing to file alphaQuartz.xyz");
    }

// supercell generation next
    try
    { 
      // Create a new file output stream
      out = new FileOutputStream("supercell.xyz");
      
      // Connect print stream to the output stream
      p = new PrintStream( out );
      
      p.println ( (num_atoms_Si + num_atoms_Ox)*(nXcell+1)*(nYcell+1)*(nZcell+1) +
                   num_atoms_C*(nXcell+1)*(nYcell+1) );    
      p.println();
// carbon
      for(iatom=0 ; iatom < num_atoms_C; iatom++ ) {
        for( i=0 ; i<=nXcell; i++ ) {
          for( j=0 ; j<=nYcell; j++ ) {
              k=nZcell;
//            Transform in cartesian coordinates
              for( icoord=0; icoord <=2; icoord++) {
                 tempo[icoord] = (carbone[iatom][0] + i)*base[0][icoord] 
                               + (carbone[iatom][1] + j)*base[1][icoord]
                               + (carbone[iatom][2] + k)*base[2][icoord];
              }
              p.printf ("C    %2.8f  %2.8f  %2.8f\n",tempo[0]*aa, tempo[1]*aa, tempo[2]*aa);
          }
        }
      }
// oxygen
      for(iatom=0 ; iatom < num_atoms_Ox; iatom++ ) {
        for( i=0 ; i<=nXcell; i++ ) {
          for( j=0 ; j<=nYcell; j++ ) {
            for( k=0 ; k<=nZcell; k++ ) {
//            Transform in cartesian coordinates
              for( icoord=0; icoord <=2; icoord++) {
                 tempo[icoord] = (oxygene[iatom][0] + i)*base[0][icoord] 
                               + (oxygene[iatom][1] + j)*base[1][icoord]
                               + (oxygene[iatom][2] + k)*base[2][icoord];
              }
              p.printf ("O    %2.8f  %2.8f  %2.8f\n",tempo[0]*aa, tempo[1]*aa, tempo[2]*aa);
            }
          }
        }
      }
// silicon
      for(iatom=0 ; iatom < num_atoms_Si; iatom++ ) {
        for( i=0 ; i<=nXcell; i++ ) {
          for( j=0 ; j<=nYcell; j++ ) {
            for( k=0 ; k<=nZcell; k++ ) {
//            Transform in cartesian coordinates
              for( icoord=0; icoord <=2; icoord++) {
                 tempo[icoord] = (silicon[iatom][0] + i)*base[0][icoord] 
                               + (silicon[iatom][1] + j)*base[1][icoord]
                               + (silicon[iatom][2] + k)*base[2][icoord];
              }
              p.printf ("Si   %2.8f  %2.8f  %2.8f\n",tempo[0]*aa, tempo[1]*aa, tempo[2]*aa);
            }
          }
        }
      }
      p.close();
    }
    catch (Exception e)
    {
      System.err.println ("Error writing to file");
    }

  }
}