#include <iostream>
#include <string.h>
#include <fstream>
#include <iomanip>
#include <sstream>
#include <stdlib.h>
#include <stdio.h>
#include <sys/stat.h>
#include <regex.hpp>
#include <regex/config.hpp>
#include <algorithm/string/split.hpp>
#include <algorithm/string.hpp>
#include <algorithm/string/regex.hpp>
#include <lexical_cast.hpp>
#include <cerrno>

/*
This test computes the unrelaxed and relaxed formation energies of a monovacancy in diamond silicon by
using the Polak-Ribiere conjugate gradient method provided in the LAMMPS 'minimize' function (lammps.sandia.gov/doc/minimize.html).
The atomic positions and the forces acting on the atoms in the unrelaxed and relaxed configurations are reported
along with the unrelaxed and relaxed formation energies.

Author: Daniel S. Karls, University of Minnesota (karl0100 |AT| umn DOT edu)
*/

using namespace std;
typedef vector<string> split_vector_type;

string get_file_contents(const char *filename){
    ifstream in(filename, std::ios::in | std::ios::binary);
    if(in){
        string contents;
        in.seekg(0, ios::end);
        contents.resize(in.tellg());
        in.seekg(0, ios::beg);
        in.read(&contents[0], contents.size());
        in.close();
        return(contents);
    }
    throw(errno);
}

int main(){
ofstream writeLAMMPS, writeEDN;
stringstream EDN_results;
int num_tot_atoms, atoms_in_unit_cell=8, tot_atom_count;
int num_cells_x, num_cells_y, num_cells_z;
double dia_vec1[3]={0.0}, dia_vec2[3]={0.0}, dia_vec3[3]={0.0}, dia_atom[8][3];
double lattice_const, cohesive_energy, xlo, xhi, ylo, yhi, zlo, zhi;
double relaxed_poteng, unrelaxed_poteng, final_pressure;
double* atom_vec;
char model[256], path[256], dblbuffer[256], dumc;
// Regular expressions
const char *relaxed_poteng_pattern="([0-9e\\.-]+)\\s+([0-9e\\.-]+)\\s\\nLoop time of ";
const char *unrelaxed_poteng_pattern="PotEng Press \\n\\s+([0-9e\\.-]+)\\s+([0-9e\\.-]+)\\s\\n";
// Use the following regex to find the relaxed atomic positions and forces.  Minimized nesting and disallowed backtracking (otherwise catastrophic
// backtracking will lead boost to run out of stack space while calling regex_search.
const char *unrelaxed_atoms_pattern="([0-9e\\.-]+\\s[0-9e\\.-]+\\n){3}(ITEM: ATOMS id type x y z fx fy fz \\n)(.*?)(ITEM: TIMESTEP)";
const char *relaxed_atoms_pattern="([0-9e\\.-]+\\s[0-9e\\.-]+\\n){3}(ITEM: ATOMS id type x y z fx fy fz \\n)(?>\\d+\\s+\\d+\\s+[0-9e\\.-]+\\s+[0-9e\\.-]+\\s+[0-9e\\.-]+\\s+[0-9e\\.-]+\\s+[0-9e\\.-]+\\s+[0-9e\\.-]+\\s+\\n)+\\z";
string output_log, dumpfile, relaxed_potengstr, unrelaxed_potengstr, relaxed_pos_and_forces, unrelaxed_pos_and_forces, stringtest;
split_vector_type split_unrelaxed_poteng, split_relaxed_poteng, split_relaxed_pos_and_forces, split_unrelaxed_pos_and_forces, split_pos;
split_vector_type split_A, split_B, split_C, split_a, split_b, split_c;
boost::cmatch relaxed_poteng_match, unrelaxed_poteng_match, relaxed_match, unrelaxed_match;
boost::regex regexpoteng_relaxed(relaxed_poteng_pattern), regexpoteng_unrelaxed(unrelaxed_poteng_pattern), regexdump_relaxed(relaxed_atoms_pattern), regexdump_unrelaxed(unrelaxed_atoms_pattern);
string atdata;

// Read from pipeline.stdin.tpl
cout<<"Please enter a valid KIM model name:"<<endl;
cin>>model;
cout<<"Please enter the equilibrium lattice constant (Angstroms) and the cohesive energy (eV) in the form [a,E_coh] (including the brackets and the comma):"<<endl;
cin>>dumc>>lattice_const>>dumc>>cohesive_energy;

// Error checks
if(lattice_const <= 0.0){
cout<<"Error: the lattice constant inputted is not positive.  Exiting..."<<endl;
cerr<<"Error: the lattice constant inputted is not positive. Exiting..."<<endl;
exit(1);
}

// Build a diamond lattice with the atom at (0,0,0) as the vacancy
dia_vec1[0]=lattice_const; dia_vec1[1]=0.0;           dia_vec1[2]=0.0;
dia_vec2[0]=0.0;           dia_vec2[1]=lattice_const; dia_vec2[2]=0.0;
dia_vec3[0]=0.0;           dia_vec3[1]=0.0;           dia_vec3[2]=lattice_const;

dia_atom[0][0]=0.0;                 dia_atom[0][1]=0.0;                 dia_atom[0][2]=0.0;
dia_atom[1][0]=0.25*lattice_const;  dia_atom[1][1]=0.25*lattice_const;  dia_atom[1][2]=0.25*lattice_const;
dia_atom[2][0]=0.0;                 dia_atom[2][1]=0.5*lattice_const;   dia_atom[2][2]=0.5*lattice_const;
dia_atom[3][0]=0.25*lattice_const;  dia_atom[3][1]=0.75*lattice_const;  dia_atom[3][2]=0.75*lattice_const;
dia_atom[4][0]=0.5*lattice_const;   dia_atom[4][1]=0.0;                 dia_atom[4][2]=0.5*lattice_const;
dia_atom[5][0]=0.75*lattice_const;  dia_atom[5][1]=0.25*lattice_const;  dia_atom[5][2]=0.75*lattice_const;
dia_atom[6][0]=0.5*lattice_const;   dia_atom[6][1]=0.5*lattice_const;   dia_atom[6][2]=0.0;
dia_atom[7][0]=0.75*lattice_const;  dia_atom[7][1]=0.75*lattice_const;  dia_atom[7][2]=0.25*lattice_const;

num_cells_x=2;
num_cells_y=num_cells_x;
num_cells_z=num_cells_x;
num_tot_atoms=num_cells_x*num_cells_y*num_cells_z*atoms_in_unit_cell-1;
xlo=0.0;
xhi=num_cells_x*lattice_const;
ylo=0.0;
yhi=num_cells_y*lattice_const;
zlo=0.0;
zhi=num_cells_z*lattice_const;
atom_vec=new double[3*num_tot_atoms];
tot_atom_count=0;

for(int xcount=0; xcount<num_cells_x; ++xcount){
    for(int ycount=0; ycount<num_cells_y; ++ycount){
        for(int zcount=0; zcount<num_cells_z;++zcount){
             for(int atomcount=0; atomcount<atoms_in_unit_cell;++atomcount){
                 if(!(xcount==0 && ycount==0 && zcount==0 && atomcount==0)){
                     atom_vec[3*tot_atom_count+0]=dia_atom[atomcount][0]+xcount*dia_vec1[0]+ycount*dia_vec2[0]+zcount*dia_vec3[0];
                     atom_vec[3*tot_atom_count+1]=dia_atom[atomcount][1]+xcount*dia_vec1[1]+ycount*dia_vec2[1]+zcount*dia_vec3[1];
                     atom_vec[3*tot_atom_count+2]=dia_atom[atomcount][2]+xcount*dia_vec1[2]+ycount*dia_vec2[2]+zcount*dia_vec3[2];
                     tot_atom_count++;
                 }
             }
        }
    }
}

// Create subdirectories under Test Result directory
mkdir("output/lammps_inputs/", S_IRWXU | S_IRWXG | S_IROTH | S_IXOTH);
mkdir("output/lammps_output_log/", S_IRWXU | S_IRWXG | S_IROTH | S_IXOTH);
mkdir("output/lammps_dump/", S_IRWXU | S_IRWXG | S_IROTH | S_IXOTH);

// Create LAMMPS input file
writeLAMMPS.open("output/lammps_inputs/in.lammps", ios::out);
if(!writeLAMMPS.fail()){
    writeLAMMPS<<"units metal\natom_style atomic\nboundary p p p\nbox tilt large\nread_data output/lammps_inputs/lammps.data"<<endl;
    writeLAMMPS<<"pair_style kim KIMvirial "<<model<<"\npair_coeff * * Si\ndump dumpid all custom 1 output/lammps_dump/lammps.dump id type x y z fx fy fz"<<endl;
    writeLAMMPS<<"dump_modify dumpid sort id format line \"\%d \%d \%16.7f \%16.7f \%16.7f \%16.7f \%16.7f \%16.7f\""<<endl;
    writeLAMMPS<<"fix 1 all box/relax iso 0.0"<<endl;
    writeLAMMPS<<"thermo_style custom pe press"<<endl;
    writeLAMMPS<<"minimize 0.00000001 0.00000001 2000 100000";
}
writeLAMMPS.close();

// Create LAMMPS data file
writeLAMMPS.open("output/lammps_inputs/lammps.data", ios::out);
if(!writeLAMMPS.fail()){
    writeLAMMPS<<"Relaxation of monovacancy in diamond silicon"<<endl<<endl;
    writeLAMMPS<<num_tot_atoms<<" atoms"<<endl<<xlo<<" "<<xhi<<" xlo xhi"<<endl<<ylo<<" "<<yhi<<" ylo yhi"<<endl<<zlo<<" "<<zhi<<" zlo zhi"<<endl<<endl;
    writeLAMMPS<<"1 atom types"<<endl<<endl<<"Atoms"<<endl<<endl;
    for(int i=0;i<num_tot_atoms;++i){
        writeLAMMPS<<i+1<<" 1 "<<atom_vec[3*i+0]<<" "<<atom_vec[3*i+1]<<" "<<atom_vec[3*i+2]<<endl;
    }
    writeLAMMPS<<endl<<"Masses"<<endl<<endl<<"1 28.0855";
}
writeLAMMPS.close();

// Free dynamic memory
delete [] atom_vec;

// Run lammps
system("lammps -in output/lammps_inputs/in.lammps > output/lammps_output_log/lammps.log");

// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// ~~~~~~~~~~  Parse lammps output & lammps dumpfile with regex ~~~~~~~~~~~
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// ******* First, parse the lammps output log for the potential energy of the unrelaxed and relaxed configurations *******
output_log=get_file_contents("output/lammps_output_log/lammps.log");
boost::regex_search(output_log.c_str(),unrelaxed_poteng_match,regexpoteng_unrelaxed);
unrelaxed_poteng=boost::lexical_cast<double>(unrelaxed_poteng_match[1]);
if(unrelaxed_poteng >= 0.0){
cout<<"Error: the unrelaxed potential energy computed in LAMMPS is not negative.  This indicates an instability of some sort has occurred.  Exiting..."<<endl;
cerr<<"Error: the unrelaxed potential energy computed in LAMMPS is not negative.  This indicates an instability of some sort has occurred.  Exiting..."<<endl;
exit(1);
}

boost::regex_search(output_log.c_str(),relaxed_poteng_match,regexpoteng_relaxed);
relaxed_potengstr = relaxed_poteng_match.str();
split(split_relaxed_poteng, relaxed_potengstr, boost::is_any_of("\n "), boost::token_compress_on);
relaxed_poteng = boost::lexical_cast<double>(split_relaxed_poteng[0]);
final_pressure = boost::lexical_cast<double>(split_relaxed_poteng[1]);
if(relaxed_poteng >= 0.0){
cout<<"Error: the relaxed potential energy computed in LAMMPS is not negative.  This indicates that an instability of some sort has occurred.  Exiting..."<<endl;
cerr<<"Error: the relaxed potential energy computed in LAMMPS is not negative.  This indicates that an instability of some sort has occurred.  Exiting..."<<endl;
exit(1);
}

// Begin writing EDN dictionary of results
EDN_results<<"{\n  \"property-id\" \"tag:staff@noreply.openkim.org,2014-04-15:property/monovacancy-formation-energy-monoatomic-cubic-diamond\"\n"
<<"  \"instance-id\" 1\n"
<<"  \"species\" {\n    \"source-value\"  \"Si\"\n  }\n"
<<"  \"a\" {\n    \"source-value\"  "<<lattice_const<<"\n    \"source-unit\"  \"angstrom\"\n  }\n"
<<"  \"unrelaxed-configuration-positions\" {\n    \"source-value\"  [\n";
// ******* Next, parse the lammps dump file for the box dimensions, atomic positions, and forces for both the relaxed and unrelaxed configurations *******
dumpfile=get_file_contents("output/lammps_dump/lammps.dump");

// Grab the unrelaxed box dimensions, atomic positions, and atomic forces.
boost::regex_search(dumpfile.c_str(),unrelaxed_match,regexdump_unrelaxed);
unrelaxed_pos_and_forces = unrelaxed_match.str();
split(split_unrelaxed_pos_and_forces, unrelaxed_pos_and_forces, boost::is_any_of("\n"), boost::token_compress_on);

// Concatenate unrelaxed positions onto EDN output file.
for(int i=0;i<num_tot_atoms;++i){
    split(split_pos, split_unrelaxed_pos_and_forces[4+i], boost::is_any_of(" "), boost::token_compress_on);
    EDN_results<<"      ["<<setw(12)<<right<<split_pos[2]<<setw(20)<<right<<split_pos[3]<<setw(20)<<right<<split_pos[4]<<"  ]\n";
    split_pos.clear();
}
EDN_results<<"    ]\n    \"source-unit\"  \"angstrom\"\n  }"<<endl;
EDN_results<<"  \"unrelaxed-configuration-forces\" {\n    \"source-value\"  [\n";
// Concatenate unrelaxed forces onto EDN output file
for(int i=0;i<num_tot_atoms;++i){
    split(split_pos, split_unrelaxed_pos_and_forces[4+i], boost::is_any_of(" "), boost::token_compress_on);
    EDN_results<<"      ["<<setw(12)<<right<<split_pos[5]<<setw(20)<<right<<split_pos[6]<<setw(20)<<right<<split_pos[7]<<" ]\n";
    split_pos.clear();
}
EDN_results<<"    ]\n    \"source-unit\"  \"eV/angstrom\"\n  }"<<endl;
EDN_results<<"  \"unrelaxed-formation-energy\" {\n"
<<"    \"source-value\"  "<<(unrelaxed_poteng+num_tot_atoms*cohesive_energy)<<"\n    \"source-unit\"  \"eV\"\n  }\n"
<<"  \"relaxed-configuration-positions\" {\n    \"source-value\"  [\n";

// Now grab the relaxed box dimensions, atomic positions, and atomic forces
boost::regex_search(dumpfile.c_str(),relaxed_match,regexdump_relaxed,boost::match_not_dot_newline);
relaxed_pos_and_forces = relaxed_match.str();
split(split_relaxed_pos_and_forces, relaxed_pos_and_forces, boost::is_any_of("\n"), boost::token_compress_on);

// Concatenate relaxed positions onto EDN output file.
for(int i=0;i<num_tot_atoms;++i){
    split(split_pos, split_relaxed_pos_and_forces[4+i], boost::is_any_of(" "), boost::token_compress_on);
    EDN_results<<"      ["<<setw(12)<<right<<split_pos[2]<<setw(20)<<right<<split_pos[3]<<setw(20)<<right<<split_pos[4]<<" ]\n";
    split_pos.clear();
}
EDN_results<<"    ]\n    \"source-unit\"  \"angstrom\"\n  }"<<endl;
EDN_results<<"  \"relaxed-configuration-forces\" {\n    \"source-value\"  [\n";

// Concatenate relaxed forces onto EDN output file
for(int i=0;i<num_tot_atoms;++i){
    split(split_pos, split_relaxed_pos_and_forces[4+i], boost::is_any_of(" "), boost::token_compress_on);
    EDN_results<<"      ["<<setw(12)<<right<<split_pos[5]<<setw(20)<<right<<split_pos[6]<<setw(20)<<right<<split_pos[7]<<" ]\n";
    split_pos.clear();
}
EDN_results<<"    ]\n    \"source-unit\"  \"eV/angstrom\"\n  }"<<endl;
EDN_results<<"  \"relaxed-formation-energy\" {\n"
<<"    \"source-value\"  "<<(relaxed_poteng+num_tot_atoms*cohesive_energy)<<"\n    \"source-unit\"  \"eV\"\n  }\n}";

// Write results to file
writeEDN.open("output/results.edn",ios::out);
writeEDN<<EDN_results.str();
writeEDN.close();

// Empty string stream
EDN_results.clear(); // Clear any set bits
EDN_results.str(string());

return 0;
}