This is a very strange problem... removing the cout in the function below causes it to stop printing the correct/expected results and printing garbage values. (i.e. it still RUNS the data it outputs, though, is wrong). Any ideas?
bool extract_tension(std::vector<double> &interfacial_tension_trap,
std::vector<double> &interfacial_tension_simp,
const std::string data,
const unsigned int num_slabs,
const double z_min, const double z_max)
{
//start @ first number
unsigned int start = 17;
unsigned int end = 17;
std::vector<double> px;
std::vector<double> py;
std::vector<double> pz;
std::vector<double> pn_minus_pt;
double result_simp=0.0;
double result_trap=0.0;
//skip timestep entry
end=get_next_space(start, data);
for(unsigned int counter=0; counter<num_slabs;counter++)
{
start = end+2;
end=get_next_space(start, data);
px.push_back(atof(data.substr(start,(end-start+1)).c_str()));
//skip the space
start = end+2;
end=get_next_space(start, data);
py.push_back(atof(data.substr(start,(end-start+1)).c_str()));
//skip the space
start = end+2;
end=get_next_space(start, data);
pz.push_back(atof(data.substr(start,(end-start+1)).c_str()));
//calculate pressure difference
// WARNING : Unit conversion ahead
// NAMD outputs pressure in bars and distance in Angstroms
// we want an integrated result of mN/m, instead.
// 1 Angstrom = 1e-10 m
// 1 bar = 1e8 mN/m^2
// net conversion -- 1e-2
pn_minus_pt.push_back((pz[counter]-0.5*(px[counter]+py[counter]))*0.01);
std::cout << "Current del_P : "
<< (pz[counter]-0.5*(px[counter]+py[counter]))*0.01
<< std::endl;
}
calculate_trapezoid(pn_minus_pt, num_slabs, z_min, z_max, result_trap);
interfacial_tension_trap.push_back(result_trap);
calculate_simpson(pn_minus_pt, num_slabs, z_min, z_max, result_simp);
interfacial_tension_simp.push_back(result_simp);
}
Apparently just touching any of vectors with a print statement allows the program to execute correctly (i.e. a printout involving px, py, OR pz)
Here's the full program:
/*********************************
*
* NAME: Interfacial Tension Calculator
* VERSION: 0.1
* AUTHOR: Jason R. Mick
* GROUP: Wayne State University, Potoff Group
* COPYRIGHT: (c) Jason R. Mick 2010
* DATE: August 9, 2010
*
* CHANGE LOG
* VERSION DATE COMMENTS
*----------------------------------------------------------------------
* 0.1 Aug. 9, 2010 Finished basic code, sans debugging
* 0.5 Aug 10, 2010 Compiled and tested code fixed error in Simpson's
* method where results were being divided rather
* than multiplied.
*
*
* FULL NOTES:
*----------------------------------------------------------------------
* You can compile this program by typing:
* g++ main.cc -o it_util
*
* You can run this program by typing:
* it_util <filename>.log <# slabs> <z-min> <z-max>
*
* where z-min and z-max represent the z-axis boundaries of the system,
* e.g.--
* it_util my_file.log 140 0.0 80.0
*
* This program only works with NAMD *.log file output
* The pressure profile MUST be turned on in your *.conf file
* for the pressure profile info to dump to the *.log file. This
* program requires that info.
*
* This program can handle 1,000+ slabs, but it has a limit to the
* character buffer and thus VERY large slab counts may cause it to fail.
*
* A standard Composite Simpson numerical integration method is used,
* which assumes a non-smooth data set.
*
* The interfacial tension is integrated at each step and then averaged
* so pertinent statistics can be gathered.
*
* You can redirect the output to store the interfacial tension
* statistics as follows:
* it_util <filename>.log <# slabs> <z-min> <z-max> > <my_file>.out
*
*******************************************/
#include <stdio.h>
#include <math.h>
#include <iostream>
#include <vector>
#include <fstream>
#include <sys/stat.h>
//Turn on to enable all interfacial
//tension results to be printed, pre-averaging
//#define DEBUG true
void start_integrations(const std::string filename,
const unsigned int num_slabs,
const double z_min, const double z_max);
int main ( int argc, char *argv[] )
{
struct stat file_info;
std::string filename = argv[1];
int slab_count;
double z_min;
double z_max;
if ( argc != 5 ) /* argc should be 3 for correct execution */
{
/*Print out proper args syntax */
std::cout << "ERROR: Missing arguments!" << std::endl
<< "Proper syntax:" << std::endl
<< "it_util <my_file>.log <# of slabs> <z-coord start>"
<< "<z-coord end>"
<< std::endl;
}
if(stat(argv[1],&file_info)==0)
{
try
{
slab_count = atoi(argv[2]);
if (slab_count > 2)
{
try
{
z_min = atof(argv[3]);
try
{
z_max = atof(argv[4]);
start_integrations(filename,
static_cast<unsigned int>(slab_count),
z_min,
z_max);
}
catch( char * str )
{
/*invalid integer third input*/
std::cout << "Invalid input -- fourth argument was invalid "
<< "decimal number, should be standard " << std::endl
<< "decimal type entry..." << std::endl
<< "I.E." << std::endl
<< "it_util my_file.log 140 0.0 80.0" << std::endl;
}
}
catch( char * str )
{
/*invalid integer third input*/
std::cout << "Invalid input -- third argument was invalid "
<< "decimal number, should be standard " << std::endl
<< "decimal type entry..." << std::endl
<< "I.E." << std::endl
<< "it_util my_file.log 140 0.0 80.0" << std::endl;
}
}
else
{
/*invalid integer secondary input*/
std::cout << "Invalid input -- second argument was invalid integer, "
<< "should be unsigned integer 2 or greater..." << std::endl
<< "I.E." << std::endl
<< "it_util my_file.log 140 0.0 80.0" << std::endl;
}
}
catch( char * str )
{
/*non integer secondary input*/
std::cout << "Invalid input -- second argument was non-integer, "
<< "should be unsigned integer 2 or greater..." << std::endl
<< "I.E." << std::endl
<< "it_util my_file.log 140 0.0 80.0" << std::endl;
}
}
else
{
/*invalid filename case...*/
std::cout << "File " << filename << "does not exist!" << std::endl
<< "Please choose valid file!" << std::endl;
}
return 1;
}
bool calculate_simpson(const std::vector<double> my_values,
const unsigned int num_points,
const double x_min, const double x_max,
double &results)
{
bool ret_val = false;
bool is_even = true;
double h;
if (my_values.size() >= 2)
{
h = (x_max-x_min)/num_points;
results+=my_values.front();
for (unsigned int counter=1; counter<num_points-1;counter++)
{
if (is_even)
{
results+=4*my_values[counter];
}
else
{
results+=2*my_values[counter];
}
is_even = !is_even;
}
results+=my_values.back();
results*=(h/3);
ret_val=true;
}
return ret_val;
}
bool calculate_trapezoid(const std::vector<double> my_values,
const unsigned int num_points,
const double x_min, const double x_max,
double &results)
{
bool ret_val = false;
double x_incr = (x_max-x_min)/(num_points-1);
if (my_values.size() >= 2)
{
for (unsigned int counter=1; counter<num_points-1; counter++)
{
results+=(x_incr/2)*(my_values[counter]+my_values[counter-1]);
}
}
return ret_val;
}
unsigned int get_next_space(const unsigned int start,
const std::string data)
{
unsigned int counter=start;
while (data.length() > counter &&
data.substr(counter,1).compare(" ") != 0)
{
counter++;
}
//if end of string, add one
if ( data.length() == counter)
counter++;
return (counter-1);
}
bool extract_tension(std::vector<double> &interfacial_tension_trap,
std::vector<double> &interfacial_tension_simp,
const std::string data,
const unsigned int num_slabs,
const double z_min, const double z_max)
{
//start @ first number
unsigned int start = 17;
unsigned int end = 17;
std::vector<double> px;
std::vector<double> py;
std::vector<double> pz;
std::vector<double> pn_minus_pt;
double result_simp=0.0;
double result_trap=0.0;
//skip timestep entry
end=get_next_space(start, data);
for(unsigned int counter=0; counter<num_slabs;counter++)
{
start = end+2;
end=get_next_space(start, data);
px.push_back(atof(data.substr(start,(end-start+1)).c_str()));
//skip the space
start = end+2;
end=get_next_space(start, data);
py.push_back(atof(data.substr(start,(end-start+1)).c_str()));
//skip the space
start = end+2;
end=get_next_space(start, data);
pz.push_back(atof(data.substr(start,(end-start+1)).c_str()));
//calculate pressure difference
// WARNING : Unit conversion ahead
// NAMD outputs pressure in bars and distance in Angstroms
// we want an integrated result of mN/m, instead.
// 1 Angstrom = 1e-10 m
// 1 bar = 1e8 mN/m^2
// net conversion -- 1e-2
pn_minus_pt.push_back((pz[counter]-0.5*(px[counter]+py[counter]))*0.01);
std::cout << "Current del_P : "
<< (pz[counter]-0.5*(px[counter]+py[counter]))*0.01
<< std::endl;
}
calculate_trapezoid(pn_minus_pt, num_slabs, z_min, z_max, result_trap);
interfacial_tension_trap.push_back(result_trap);
calculate_simpson(pn_minus_pt, num_slabs, z_min, z_max, result_simp);
interfacial_tension_simp.push_back(result_simp);
}
double average_vector(std::vector<double> my_vector)
{
double average_val=0.0;
for(unsigned int counter=0; counter< my_vector.size(); counter++)
{
average_val+=my_vector[counter]/my_vector.size();
}
return average_val;
}
double std_dev_vector(std::vector<double> my_vector)
{
double std_deviation=0.0;
double average_val = average_vector(my_vector);
for(unsigned int counter=0; counter< my_vector.size(); counter++)
{
std_deviation+=(my_vector[counter]-average_val)*
(my_vector[counter]-average_val);
}
std_deviation=sqrt(std_deviation);
return std_deviation;
}
void start_integrations(const std::string filename,
const unsigned int num_slabs,
const double z_min, const double z_max)
{
std::ifstream in_file;
std::vector<double> interfacial_tension_trap;
std::vector<double> interfacial_tension_simp;
std::string current_line;
char * cstr_line;
bool data_grab_success = true;
in_file.open(filename.c_str(), std::ifstream::in);
while (!in_file.eof() && data_grab_success)
{
cstr_line=(char *) malloc(sizeof(char)*65536);
//get new line
in_file.getline(cstr_line,65536);
current_line = cstr_line;
free(cstr_line);
if (current_line.substr(0,15).compare("PRESSUREPROFILE")==0)
{
//pressure profile found!
//process line to get the interfacial tension, check that it succeeded
data_grab_success = extract_tension(interfacial_tension_trap,
interfacial_tension_simp,
current_line,
num_slabs,
z_min,
z_max);
}
}
in_file.close();
//print stats
std::cout << "Interfacial Tension (Trapezoid Method): "
<< average_vector(interfacial_tension_trap) << std::endl
<< "Standard Deviation (Trapezoid Method): "
<< std_dev_vector(interfacial_tension_trap) << std::endl
<< "Interfacial Tension (Composite Simpson's Method): "
<< average_vector(interfacial_tension_simp) << std::endl
<< "Standard Deviation (Composite Simpson's Method): "
<< std_dev_vector(interfacial_tension_simp) << std::endl;
}
And here's a sample set of data:
Removed... see explanation at end of post for link to data to use.
Compile like so:
g++ main.cc -o it_util
Run using the command:
it_util equil2_NVT_PP_318Slabs.log 318 0.0 318.0 > temp.out
FYI before someone comments on my #ifdef "debug" statements, please note that they are for data dumping. I have used GDB before. I'm guessing had I not said this, someone will comment "Learn to use gdb." In this case the program loops through so many iterations, GDB doesn't give me useful info, where printouts dumped to a output file DO.
NOTE:
Actually I discovered that if you use the chopped down version of the file being parsed (in the data section above) the program also doesn't output the correct data. When I restored the original data file it worked, but the file is too large to post here (I tried...) so that is not an option....
Instead I've uploaded a full pastebin to here:
http://pastebin.com/JasbSc7B