#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <iostream>
#include <fstream>
#include <stdexcept>
#include <sstream>
#include <algorithm>
#include <string>
#include <vector>
#include <map>

#include <iomanip> 
using namespace std;

#include <boost/bind.hpp>
#include <boost/function.hpp>

#include <boost/iostreams/device/file.hpp>
#include <boost/iostreams/filtering_stream.hpp>
#include <boost/iostreams/filter/gzip.hpp>

using namespace boost;
using namespace boost::iostreams;


 int main()
   {
	double rcut = 1.2;

	double krf = 0.289352;
	double crf = 1.25; 
	double epsfac = 9.26236;	
	
	ofstream pf("pf.log");
	
	unsigned int N = 1000;
	double delt = rcut/double(N);
	std::vector<double> potential;
	potential.resize(N);
	std::vector<double> force;
	force.resize(N);

	for (unsigned int i =1; i< N; i++)
		{
		double r = double(i)*delt;
		double rsq = r*r;
		double rinv = 1.0/r;
		double r2inv = rinv*rinv;
	
		double force_divr = epsfac * (r2inv - double(2.0)*krf*r);
		double pair_eng = epsfac * (rinv + krf*rsq - crf);			
	
		force[i] = force_divr;
		potential[i] = pair_eng;			
		} 
		
	for (unsigned int i =1; i< N; i++)
		{
		double r = double(i)*delt;
		pf<<r<<"  "<<potential[i]<<"  "<<force[i]<<endl;
		}
	}