#include"afterprocess.h"

namespace EVO {
	double get_energy()
	{
		double EDIFFG = atof(GetInfoINCAR("EDIFFG").c_str()); //EDIFFG < 0
		double energy = 0;
		POSCAR pos;
		FILE* fp = fopen("POSCAR", "r");
		if (fp == NULL)
		{
			printf("POSCAR IS NOT EXIST!\n");
			return energy;
		}
		readposcar(fp, pos);
		fclose(fp);
		FILE* fp_ = fopen("OUTCAR", "r");
		if (fp_ == NULL)
		{
			printf("OUTCAR IS NOT EXIST!\n");
			return energy;
		}
		char buf[1024];
		bool flag1 = 0, flag2 = 0;
		vector<vector<double> > force(pos.nant[0], vector<double>(3));
		while (fgets(buf, 1024, fp_) != NULL)
		{
			if (strstr(buf, "energy(sigma->0)"))
			{
				sscanf(buf, "%*s%*s%*s%*s%*s%*s%*s%lf", &energy);
			}
			if (strstr(buf, "TOTAL-FORCE"))
			{
				flag1 = 1;
				fgets(buf, 1024, fp_); // ------------------
				for (int i = 0; i < pos.nant[0]; i++)
				{
					fgets(buf, 1024, fp_);
					sscanf(buf, "%*lf%*lf%*lf%lf%lf%lf", &force[i][0], &force[i][1], &force[i][2]);
				}
			}
		}
		fclose(fp_);
		if (flag1 == 1)
		{
			for (int i = 0; i < pos.nant[0]; i++)
				for (int j = 0; j < 3; j++)
				{
					if (fabs(force[i][j]) > fabs(EDIFFG))
					{
						flag2 = 1;
					}
				}
		}
		else
			flag2 = 1;
		return flag2 == 1 ? 0 : energy;
	}

	double get_volume()
	{
		FILE* fp = fopen("CONTCAR", "r");
		if (fp == NULL)
			return 0;
		int line = 0;
		double vec[3][3];
		char buf[100];
		while (line < 5)
		{
			fgets(buf, 100, fp);
			if (line == 2)
				sscanf(buf, "%lf%lf%lf", &vec[0][0], &vec[0][1], &vec[0][2]);
			if (line == 3)
				sscanf(buf, "%lf%lf%lf", &vec[1][0], &vec[1][1], &vec[1][2]);
			if (line == 4)
				sscanf(buf, "%lf%lf%lf", &vec[2][0], &vec[2][1], &vec[2][2]);
			line++;
		}
		double odd = 0;
		double even = 0;
		for (int i = 0; i < 3; i++)
			odd += vec[0][i] * vec[1][(i + 1) % 3] * vec[2][(i + 2) % 3];
		for (int i = 0; i < 3; i++)
			even += vec[0][i] * vec[1][(i + 2) % 3] * vec[2][(i + 1) % 3];
		return fabs(odd - even);
	}
}
