/* We illustrate the operation of the code to optimize 130 atoms system
atomic coordinations are in file Ni.pdb
For default parameters, the results are (Dos environment):
==========================================
..........100  Value of evaluate function : -413.905
..........200  Value of evaluate function : -469.3
..........300  Value of evaluate function : -475.807
..........400  Value of evaluate function : -476.588
..........500  Value of evaluate function : -477.183
..........600  Value of evaluate function : -477.333
..........700  Value of evaluate function : -469.028
..........800  Value of evaluate function : -478.224
..........900  Value of evaluate function : -478.761
..........1000  Value of evaluate function : -478.735
..........1100  Value of evaluate function : -479.298
..........1200  Value of evaluate function : -479.34
..........1300  Value of evaluate function : -479.355
..........1400  Value of evaluate function : -479.447
..........1500  Value of evaluate function : -479.828
..........1600  Value of evaluate function : -479.958
..........1700  Value of evaluate function : -480.003
..........1800  Value of evaluate function : -480.016
..........1900  Value of evaluate function : -480.019
..........2000  Value of evaluate function : -480.019
..........2100  Value of evaluate function : -480.02
..........2200  Value of evaluate function : -480.019
..........2300  Value of evaluate function : -480.019
..........2400  Value of evaluate function : -480.019

Final convergence status = 0
Convergence tolerance for gradient satisfied
projected gradient max norm: 9.996547e-004
function value:              -4.800188e+002

Total cg  iterations:                  769
Total cg  function evaluations:       2161
Total cg  gradient evaluations:       1504
Total cbb iterations:                  251
Total cbb function evaluations:        507
Total cbb gradient evaluations:        252
------------------------------------------
Total function evaluations:           2668
Total gradient evaluations:           1756
==========================================
==========================================
*/

#include <math.h> 
#include "asa_user.h" /* needed by the program which calls asa_cg */
#include "../PackingGA/pair_sc.h"
#include "../PackingGA/pair_sc.cpp"
#include "../PackingGA/CParticles.h"
#include "../PackingGA/CParticles.cpp"
#include "../PackingGA/CRand.h"
#include "../PackingGA/CRand.cpp"
#include <fstream>
#include <iostream>
#include <sstream>
using namespace std;

int dim = 390;//390;
int num = dim/3;
CParticles *atoms;

/* prototypes for the function and gradient evaluation routines */
double EnergySC_CYL(double *xx, int n); /// Caution: n = num*3
void ForceSC_CYL(double  *g, double  *xx, int n); /// Caution: n = num*3
double ForceEnergySC_CYL (double  *g, double  *xx, INT32 n) ;

int main (int argc, char **argv)
{
	//double epsilon = 1.e-4;
	double grad_tol = 1.e-3;
	double g_rho_min = 0;
	double g_rho_max = 2.8;
	double g_phi_min = -PI;
	double g_phi_max = PI;
	double g_z_min = 0;
	double g_z_max = 36;
	double rho_scale = 0.78;
	double z_scale = 1.0;
	int center = 1;
	string infile("Ni.pdb");
	string outpdb("Ni_evolve.pdb");
	string outtxt("Ni_output.txt");
	string element("Ni");
	string parafile("param.ini");

	/// parse command to obtain parameters filename.
	for(int i = 1; i < argc; i++)
	{
		stringstream str;
		if(strcmp(argv[i], "-parafile") == 0){
			i++;
			if(i >= argc){
				cerr << argv[0] << ": -parafile needs a filename.\n";
				exit(1);
			}
			else{
				str << argv[i];
				str >> parafile;
			}
		}
	}

	/// Parsing param.ini
	ifstream fin(parafile.c_str());
	if(!fin)
	{
		cerr << "-------------------------------------------------------" << endl;
		cerr << "Warning: Cannot open " << parafile.c_str() << " for input. " << endl;
		cerr << "	 Default or command parameters will be used." << endl;
		cerr << "-------------------------------------------------------" << endl;
	}
	else
	{
		cout << "=======================================================" << endl;
		cout << "Now using parameters in file " << parafile.c_str() << endl;
		cout << "Only command line options can revise these parameters. " << endl;
		cout << "=======================================================" << endl;
		while(fin.good() && !fin.eof())
		{
			stringstream str;
			string line, tmp;
			getline(fin, line);
			if(line.find("grad_tol") != line.npos){
				str << line;
				str >> tmp >> grad_tol;
			}
			else if(line.find("g_rho_min") != line.npos){
				str << line;
				str >> tmp >> g_rho_min;
			}
			else if(line.find("g_rho_max") != line.npos){
				str << line;
				str >> tmp >> g_rho_max;
			}
			else if(line.find("g_phi_min") != line.npos){
				str << line;
				str >> tmp >> g_phi_min;
			}
			else if(line.find("g_phi_max") != line.npos){
				str << line;
				str >> tmp >> g_phi_max;
			}
			else if(line.find("g_z_min") != line.npos){
				str << line;
				str >> tmp >> g_z_min;
			}
			else if(line.find("g_z_max") != line.npos){
				str << line;
				str >> tmp >> g_z_max;
			}
			else if(line.find("rho_scale") != line.npos){
				str << line;
				str >> tmp >> rho_scale;
			}
			else if(line.find("z_scale") != line.npos){
				str << line;
				str >> tmp >> z_scale;
			}
			else if(line.find("infile") != line.npos){
				str << line;
				str >> tmp >> infile;
			}
			else if(line.find("outpdb") != line.npos){
				str << line;
				str >> tmp >> outpdb;
			}
			else if(line.find("outtxt") != line.npos){
				str << line;
				str >> tmp >> outtxt;
			}
			else if(line.find("center") != line.npos){
				str << line;
				str >> tmp >> center;
			}
			else if(line.find("element") != line.npos){
				str << line;
				str >> tmp >> element;
			}
			else{
				cout << "Parameters in file param.ini are wrong! \n";
				cout << "Default or command paramters will be used. \n";
			}
		}
	}
	fin.close();

	for(int i = 1; i < argc; i++)
	{
		stringstream str;
		if(strcmp(argv[i], "-grad_tol") == 0){
			i++;
			if(i >= argc){
				cerr << argv[0] << ": -grad_tol needs a parameter.\n";
				exit(1);
			}
			else{
				str << argv[i];
				str >> grad_tol;
				continue;
			}
		}
		else if(strcmp(argv[i], "-g_rho_min") == 0){
			i++;
			if(i >= argc){
				cerr << argv[0] << ": -g_rho_min needs a paramter.\n";
				exit(1);
			}
			else{
				str << argv[i];
				str >> g_rho_min;
				continue;
			}
		}
		else if(strcmp(argv[i], "-g_rho_max") == 0){
			i++;
			if(i >= argc){
				cerr << argv[0] << ": -g_rho_max needs a parameter.\n";
				exit(1);
			}
			else{
				str << argv[i];
				str >> g_rho_max;
				continue;
			}
		}
		else if(strcmp(argv[i], "-g_phi_min") == 0){
			i++;
			if(i >= argc){
				cerr << argv[0] << ": -g_phi_min needs a paramter.\n";
				exit(1);
			}
			else{
				str << argv[i];
				str >> g_phi_min;
			}
		}
		else if(strcmp(argv[i], "-g_phi_max") == 0){
			i++;
			if(i >= argc){
				cerr << argv[0] << ": -g_phi_max needs a paramter.\n";
				exit(1);
			}
			else{
				str << argv[i];
				str >> g_phi_max;
			}
		}
		else if(strcmp(argv[i], "-g_z_min") == 0){
			i++;
			if(i >= argc){
				cerr << argv[0] << ": -g_z_min needs a parameter.\n";
				exit(1);
			}
			else{
				str << argv[i];
				str >> g_z_min;
			}
		}
		else if(strcmp(argv[i], "-g_z_max") == 0){
			i++;
			if(i >= argc){
				cerr << argv[0] << ": -g_z_max needs a parameter.\n";
				exit(1);
			}
			else{
				str << argv[i];
				str >> g_z_max;
			}
		}
		else if(strcmp(argv[i], "-rho_scale") == 0){
			i++;
			if(i >= argc){
				cerr << argv[0] << ": -rho_scale needs a paramter.\n";
				exit(1);
			}
			else{
				str << argv[i];
				str >> rho_scale;
			}
		}
		else if(strcmp(argv[i], "-z_scale") == 0){
			i++;
			if(i >= argc){
				cerr << argv[0] << ": -z_scale needs a paramter.\n";
				exit(1);
			}
			else{
				str << argv[i];
				str >> z_scale;
			}
		}
		else if(strcmp(argv[i], "-infile") == 0){
			i++;
			if(i >= argc){
				cerr << argv[0] << ": -infile needs a filename.\n";
				exit(1);
			}
			else{
				str << argv[i];
				str >> infile;
			}