#include<fstream>
#include<iostream>
#include<sstream>
#include<string>
#include <cstdlib>
#include <iomanip>

using namespace std;

int main(int argc, char *argv[]){

if(argc != 3) {cerr << "ERROR: wrong number of arguments." << endl; return 0;}

bool OK = false;
string filename(argv[1]); // name of the dump file
ifstream infile(filename.c_str());

int ref_step = atoi(argv[2]); // time step from the dump file at the surface tention should be calculated

double Lx = 24.56;
double Ly = Lx;
double A  = 2.*Lx*Ly;	// times 2 because there is two surfaces

double m = 33.308750217; //times 10^{-26}kg, m = M/Na, M - molar mass, Na - Avogadro number

string line;
int n, time_step;
double x,y,z,vx,vy,vz,fx,fy,fz;

double gamma_vir = 0.0; // virial contribution to the pressure
double gamma_kin = 0.0; // kinetic energy contribution to the pressure

while (getline(infile,line))
{	      
  if(line.compare("ITEM: TIMESTEP") == 0) {
	getline(infile,line);
	std::istringstream iss(line);
	iss >> time_step;
	OK = false; 
	continue; 
   }
  
  if(line.compare("ITEM: ATOMS id x y z vx vy vz fx fy fz ") == 0) { OK = true; continue; }
  
  if(OK){
	if(time_step == ref_step){
		std::istringstream iss(line);
		iss >> n >> x >> y >> z >> vx >> vy >> vz >> fx >> fy >> fz;
		gamma_vir = gamma_vir + z*fz - (x*fx + y*fy)/2.;
		gamma_kin = gamma_kin + vz*vz - (vx*vx + vy*vy)/2.;
	}
  }
}

double gamma;

gamma = ( m*gamma_kin*0.01 + 16.*gamma_vir )/A/3.;  // gamma_kin is multiplied by 0.01 because of the conversion of the picoseconds into seconds in the denominator then multiplying in the numerator with 10^{-26} which comes from the mass; the second term is in the units eV/Angstroem^2 = 16 N/m
cout << "# surface tension = " << setw(12) << fixed << setprecision(9) << gamma << endl;

return 1;
}