Skip to content
Snippets Groups Projects
Commit a42b0b7d authored by Axel Kohlmeyer's avatar Axel Kohlmeyer
Browse files

update colvars abf integrate tool from git repo

parent 03828b58
No related branches found
No related tags found
No related merge requests found
...@@ -12,13 +12,12 @@ ABFdata::ABFdata(const char *gradFileName) ...@@ -12,13 +12,12 @@ ABFdata::ABFdata(const char *gradFileName)
std::ifstream gradFile; std::ifstream gradFile;
std::ifstream countFile; std::ifstream countFile;
int n;
char hash; char hash;
double xi; double xi;
char *countFileName; char *countFileName;
countFileName = new char[strlen (gradFileName) + 2]; countFileName = new char[strlen (gradFileName) + 2];
strcpy (countFileName, gradFileName); strcpy (countFileName, gradFileName);
countFileName[strlen (gradFileName) - 4] = '\0'; countFileName[strlen (gradFileName) - 4] = '\0';
strcat (countFileName, "count"); strcat (countFileName, "count");
...@@ -152,7 +151,8 @@ ABFdata::ABFdata(const char *gradFileName) ...@@ -152,7 +151,8 @@ ABFdata::ABFdata(const char *gradFileName)
} }
// Could check for end-of-file string here // Could check for end-of-file string here
countFile.close(); countFile.close();
delete [] countFileName; delete[] countFileName;
delete[] pos;
// for metadynamics // for metadynamics
bias = new double[scalar_dim]; bias = new double[scalar_dim];
...@@ -269,7 +269,7 @@ void ABFdata::write_bias(const char *fileName) ...@@ -269,7 +269,7 @@ void ABFdata::write_bias(const char *fileName)
if (minbias == 0.0 || (bias[index] > 0.0 && bias[index] < minbias)) if (minbias == 0.0 || (bias[index] > 0.0 && bias[index] < minbias))
minbias = bias[index]; minbias = bias[index];
} }
for (index = 0; index < scalar_dim; index++) { for (index = 0; index < scalar_dim; index++) {
// Here we do the Euclidian division iteratively // Here we do the Euclidian division iteratively
for (i = Nvars - 1; i > 0; i--) { for (i = Nvars - 1; i > 0; i--) {
......
...@@ -22,10 +22,7 @@ int main(int argc, char *argv[]) ...@@ -22,10 +22,7 @@ int main(int argc, char *argv[])
char *out_file; char *out_file;
unsigned int step, nsteps, total, out_freq; unsigned int step, nsteps, total, out_freq;
int *pos, *dpos, *newpos; int *pos, *dpos, *newpos;
unsigned int *histogram;
const double *grad, *newgrad;
unsigned int offset, newoffset; unsigned int offset, newoffset;
int not_accepted;
double dA; double dA;
double temp; double temp;
double mbeta; double mbeta;
...@@ -47,7 +44,7 @@ int main(int argc, char *argv[]) ...@@ -47,7 +44,7 @@ int main(int argc, char *argv[])
if (!(data_file = parse_cl(argc, argv, &nsteps, &temp, &meta, &hill, &hill_fact))) { if (!(data_file = parse_cl(argc, argv, &nsteps, &temp, &meta, &hill, &hill_fact))) {
std::cerr << "\nabf_integrate: MC-based integration of multidimensional free energy gradient\n"; std::cerr << "\nabf_integrate: MC-based integration of multidimensional free energy gradient\n";
std::cerr << "Version 20110511\n\n"; std::cerr << "Version 20160420\n\n";
std::cerr << "Syntax: " << argv[0] << std::cerr << "Syntax: " << argv[0] <<
" <filename> [-n <nsteps>] [-t <temp>] [-m [0|1] (metadynamics)]" " <filename> [-n <nsteps>] [-t <temp>] [-m [0|1] (metadynamics)]"
" [-h <hill_height>] [-f <variable_hill_factor>]\n\n"; " [-h <hill_height>] [-f <variable_hill_factor>]\n\n";
...@@ -62,7 +59,7 @@ int main(int argc, char *argv[]) ...@@ -62,7 +59,7 @@ int main(int argc, char *argv[])
} else { } else {
std::cout << "\nUsing unbiased MC sampling\n"; std::cout << "\nUsing unbiased MC sampling\n";
} }
if (nsteps) { if (nsteps) {
std::cout << "Sampling " << nsteps << " steps at temperature " << temp << "\n\n"; std::cout << "Sampling " << nsteps << " steps at temperature " << temp << "\n\n";
out_freq = nsteps / 10; out_freq = nsteps / 10;
...@@ -74,7 +71,7 @@ int main(int argc, char *argv[]) ...@@ -74,7 +71,7 @@ int main(int argc, char *argv[])
converged = false; converged = false;
} }
// Inverse temperature in (kcal/mol)-1 // Inverse temperature in (kcal/mol)-1
mbeta = -1 / (0.001987 * temp); mbeta = -1 / (0.001987 * temp);
ABFdata data(data_file); ABFdata data(data_file);
...@@ -91,12 +88,16 @@ int main(int argc, char *argv[]) ...@@ -91,12 +88,16 @@ int main(int argc, char *argv[])
dpos = new int[data.Nvars]; dpos = new int[data.Nvars];
newpos = new int[data.Nvars]; newpos = new int[data.Nvars];
// TODO: this will be an infinite loop if there is no sampling
// it would be more robust to build a list of non-empty bins, and
// pick from there (or just treat the special case of no sampling
// and output a null PMF)
do { do {
for (int i = 0; i < data.Nvars; i++) { for (int i = 0; i < data.Nvars; i++) {
pos[i] = rand() % data.sizes[i]; pos[i] = rand() % data.sizes[i];
} }
offset = data.offset(pos); offset = data.offset(pos);
} while ( !data.allowed (offset) ); } while ( !data.allowed (offset) );
rmsd = compute_deviation(&data, meta, 0.001987 * temp); rmsd = compute_deviation(&data, meta, 0.001987 * temp);
std::cout << "\nInitial gradient RMS is " << rmsd << "\n"; std::cout << "\nInitial gradient RMS is " << rmsd << "\n";
...@@ -128,9 +129,9 @@ int main(int argc, char *argv[]) ...@@ -128,9 +129,9 @@ int main(int argc, char *argv[])
data.bias[offset] += hill; data.bias[offset] += hill;
} }
grad = data.gradients + offset * data.Nvars; const double *grad = data.gradients + offset * data.Nvars;
not_accepted = 1; int not_accepted = 1;
while (not_accepted) { while (not_accepted) {
dA = 0.0; dA = 0.0;
total++; total++;
...@@ -206,7 +207,7 @@ double compute_deviation(ABFdata * data, bool meta, double kT) ...@@ -206,7 +207,7 @@ double compute_deviation(ABFdata * data, bool meta, double kT)
double *dev = data->deviation; double *dev = data->deviation;
double *est = data->estimate; double *est = data->estimate;
const double *grad = data->gradients; const double *grad = data->gradients;
int *pos, *dpos, *newpos; int *pos, *newpos;
double rmsd = 0.0; double rmsd = 0.0;
unsigned int offset, newoffset; unsigned int offset, newoffset;
double sum; double sum;
...@@ -215,7 +216,6 @@ double compute_deviation(ABFdata * data, bool meta, double kT) ...@@ -215,7 +216,6 @@ double compute_deviation(ABFdata * data, bool meta, double kT)
unsigned int norm = 0; // number of data points summmed unsigned int norm = 0; // number of data points summmed
pos = new int[data->Nvars]; pos = new int[data->Nvars];
dpos = new int[data->Nvars];
newpos = new int[data->Nvars]; newpos = new int[data->Nvars];
for (int i = 0; i < data->Nvars; i++) for (int i = 0; i < data->Nvars; i++)
...@@ -287,7 +287,6 @@ double compute_deviation(ABFdata * data, bool meta, double kT) ...@@ -287,7 +287,6 @@ double compute_deviation(ABFdata * data, bool meta, double kT)
delete [] pos; delete [] pos;
delete [] newpos; delete [] newpos;
delete [] dpos;
return sqrt(rmsd / norm); return sqrt(rmsd / norm);
} }
...@@ -296,8 +295,6 @@ double compute_deviation(ABFdata * data, bool meta, double kT) ...@@ -296,8 +295,6 @@ double compute_deviation(ABFdata * data, bool meta, double kT)
char *parse_cl(int argc, char *argv[], unsigned int *nsteps, double *temp, char *parse_cl(int argc, char *argv[], unsigned int *nsteps, double *temp,
bool * meta, double *hill, double *hill_fact) bool * meta, double *hill, double *hill_fact)
{ {
char *filename = NULL;
float f_temp, f_hill;
int meta_int; int meta_int;
// getting default value for the integer // getting default value for the integer
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment