/*Generates cluster positions for fcc and bcc 
 * noble metals Cu, Au and Ag.
 * To compile: g++ pos_gen.cpp -o gen
 * To run: ./gen 1 (for bcc clusters)
 *         ./gen 2 (for fcc clusters)
 * 
 * position files output to same directory.
 * =======
 * S. Silayi - 01/2018
 *
*/ 
#include <stdio.h>
#include <math.h>
#include <cstdlib>
#include <fstream>
#include <iostream>
#include <string>
#include <iomanip>

using namespace std;

void initPositions( int key, int N, double ** r, double * dist, int &index ) {
    
    double lc = 1.0;
    
    int nbcc = 2;
    int nfcc = 4;
    
    int nc; 
    double * dx; double *dy; double *dz;
    if (key == 1) {
        nc = nbcc;
        
        dx = new double [nc];
        dy = new double [nc];
        dz = new double [nc];
        
        // positions in fcc unit cell
        dx[0] = 0.0; dx[1] = 0.5;
        dy[0] = 0.0; dy[1] = 0.5;
        dz[0] = 0.0; dz[1] = 0.5;
        
    }
    
    else if (key == 2){
        
        nc = nfcc;   
        dx = new double [nc];
        dy = new double [nc];
        dz = new double [nc];
        
        // positions in fcc unit cell
        dx[0] = 0.0; dx[1] = 0.5;dx[2] = 0.5; dx[3] = 0.0;
        dy[0] = 0.0; dy[1] = 0.5;dy[2] = 0.0; dy[3] = 0.5;
        dz[0] = 0.0; dz[1] = 0.0;dz[2] = 0.5; dz[3] = 0.5;
        
    }
    
    // find M large enough to fit N atoms on a lattice
    int M = 0.0;
    while (nc*M*M*M < N)
        M = M+1;
    
    double L = M*lc;
    
    int n = 0;
    for (int x = 0; x < M; x++){
        for (int y = 0; y < M; y++){
            for (int z = 0; z < M; z++)
                for (int p = 0; p < nc; p++){
                    if (n < N){
                        r[n][0] = (x + dx[p]) * lc;
                        r[n][1] = (y + dy[p]) * lc;
                        r[n][2] = (z + dz[p]) * lc;
                        ++n;
                    }
                }
        }
    }
    
    
    double rCM[3] = {0, 0, 0};
    for (int n = 0; n < N; n++)
        for (int i = 0; i < 3; i++)
            rCM[i] += r[n][i];
        for (int i = 0; i < 3; i++)
            rCM[i] /= N;
        for (int n = 0; n < N; n++)
            for (int i = 0; i < 3; i++)
                r[n][i] -= rCM[i];
            
            
            double center[3] = {0.0, 0.0, 0.0};
        
        for (int i = 0; i < N; i++){
            
            
            double rSqd;
            double dr[3] = {0.0};
            
            
            dr[0] = -r[i][0] + center[0];
            dr[1] = -r[i][1] + center[1];
            dr[2] = -r[i][2] + center[2];
            
            
            rSqd = dr[0]*dr[0] + dr[1]*dr[1] + dr[2]*dr[2];
            dist[i] = sqrt(rSqd);
            
        }
        
        int min = dist[0];
        index = 0;
        for (int i = 1; i < N; i++)
            if (dist[i]<min) {min = dist[i]; index = i;}
            
            double shift_x = r[index][0];
        double shift_y = r[index][1];
        double shift_z = r[index][2];
        
        for (int i = 0; i < N; i++){
            
            r[i][0] =  (r[i][0]-shift_x);
            r[i][1] =  (r[i][1]-shift_y);
            r[i][2] =  (r[i][2]-shift_z);
            
            double dx = center[0] - r[i][0];
            double dy = center[1] - r[i][1];
            double dz = center[2] - r[i][2];
            double rSqd = dx*dx + dy*dy + dz*dz;
            
            dist[i] = sqrt(rSqd);
            
        }
        
        return;
}


void write_pos(double ** pos, int ncount, string out, double a){
    
    
    ofstream outfile; outfile.open(out.c_str());
    
    char K;
    string K1;
    int k = 0;
    int prec;
    
    while (k < ncount)
    {
        
        if (fabs(pos[k][0]*a) > 9.9) prec =3;
        else prec = 4;
        
        if (pos[k][0] <0)
            outfile<<"      "<< fixed<<setprecision(prec)<<pos[k][0]*a;
        else
            outfile<<"       "<< fixed<< setprecision(prec)<<pos[k][0]*a;
        
        if (fabs(pos[k][1]*a) > 9.9) prec =3;
        else prec = 4;
        
        if (pos[k][1] <0)
            outfile<<"    "<< fixed<< setprecision(prec)<<pos[k][1]*a;
        else 
            outfile<<"     "<< fixed<< setprecision(prec)<<pos[k][1]*a;
        
        if (fabs(pos[k][2]*a) > 9.9) prec =3;
        else prec = 4;
        
        if (pos[k][2] <0)
            outfile<<"    "<< fixed<< setprecision(prec)<<pos[k][2]*a;
        else
            outfile<<"     "<< fixed<< setprecision(prec)<<pos[k][2]*a;
        
        outfile << "\n";
        
        k=k+1;
    }
    
    
    outfile.close();
    
}




int main(int argc, char *argv[]){
    
    
    int key;
    if ( argc != 2 ) {
        cout<<"usage: "<< argv[0] <<" <key (1=bcc, 2=fcc)>\n";
        cout << "Exiting! \n";
        return 0;
    }
    else {
        
        key = atoi(argv[1]);
        if (key == 1){
        
            cout<<"Generating BCC Cluster Positions: \n";
            cout << "Cluster \t Count  \n";
        }
        else if (key ==2){
            
            cout<<"Generating FCC Cluster Positions: \n";
            cout << "Cluster \t Count  \n";
        }
    }  
    int N = 100000;
    double ** r = new double* [N];
    double ** nnr = new double* [N]; 
    double * dist = new double [N];
    int index;
    
    for (int i = 0; i < N; i++){ 
        r[i] = new double[3];
        nnr[i] = new double[3];
    }
    
    initPositions( key, N, r, dist, index );
    
    double * nn = new double [14];
    
    if (key ==1){ 
        double nnbcc[14] = {0.0, sqrt(3.0)/2.0, 1.0, sqrt(2.0), sqrt(11.0/4.0), 
            sqrt(3.0), 2.0, sqrt(19.0/4.0), sqrt(5.0), sqrt(6.0), 
            sqrt(27.0/4.0), sqrt(8.0), sqrt(35.0/4.0), 3.0}; 
            nn= nnbcc;
    }else if (key ==2){                        
        double nnfcc[14] = {0.0, sqrt(1.0/2.0), 1.0, sqrt(6.0/4.0), sqrt(2.0), 
            sqrt(10.0/4.0), sqrt(3.0), sqrt(14.0/4.0), 2.0, sqrt(18.0/4.0), 
            sqrt(5.0), sqrt(22.0/4.0), sqrt(6.0), sqrt(26.0/4.0)};                          
            
            nn = nnfcc;
    }
    int ncount = 0;
    int n_cum = 0; 
    for (int j = 0; j < 14; j++){    
        int nc = 0;
        
        for (int i =0; i < N; i++ )
        {
            
            if (dist[i] == nn[j]){
                
                nnr[ncount][0] = r[i][0];            
                nnr[ncount][1] = r[i][1];            
                nnr[ncount][2] = r[i][2];            
                
                ncount++;        
                nc++;
            }
        }

	n_cum +=nc;
        cout << j << "\t" << nc << "\t" << n_cum << "\n";
    }
    
    if (key == 1){
        
        double a_cu = 2.80;
        double a_au = 3.23; 
        double a_ag = 3.23;
        
        
        string out = "cu_bcc.pos";
        double a = a_cu;
        write_pos(nnr, ncount, out, a);
        
        out = "au_bcc.pos";
        a = a_au;
        write_pos(nnr, ncount, out, a);
        
        out = "ag_bcc.pos";
        a = a_ag;
        write_pos(nnr, ncount, out, a);
        
        a = 1.;
        out = "bcc_clusters";
        write_pos(nnr, ncount, out, a);
        cout << "Total: " << ncount << "\n";   
    }
    
    else if (key == 2){
        
        double a_cu = 3.53;
        double a_au = 4.06; 
        double a_ag = 4.01;
        
        
        string out = "cu_fcc.pos";
        double a = a_cu;
        write_pos(nnr, ncount, out, a);
        
        out = "au_fcc.pos";
        a = a_au;
        write_pos(nnr, ncount, out, a);
        
        out = "ag_fcc.pos";
        a = a_ag;
        write_pos(nnr, ncount, out, a);
        
        a = 1.;
        out = "fcc_clusters";
        write_pos(nnr, ncount, out, a);
        cout << "Total: " << ncount << "\n";   
    }
    
    return 0;    
    
}