mirror of https://github.com/lammps/lammps.git
344 lines
11 KiB
C++
344 lines
11 KiB
C++
|
/****************************************************************
|
||
|
|
* abf_integrate *
|
||
|
|
* Integrate n-dimensional PMF from discrete gradient grid *
|
||
|
|
* Jerome Henin <jerome.henin@ibpc.fr> *
|
||
|
|
****************************************************************/
|
||
|
|
|
||
|
|
#include "abf_data.h"
|
||
|
|
#include <fstream>
|
||
|
|
#include <string>
|
||
|
|
#include <cstring>
|
||
|
|
#include <cstdlib>
|
||
|
|
#include <ctime>
|
||
|
|
#include <cmath>
|
||
|
|
|
||
|
|
char *parse_cl(int argc, char *argv[], unsigned int *nsteps, double *temp,
|
||
|
|
bool * meta, double *hill, double *hill_fact);
|
||
|
|
double compute_deviation(ABFdata * data, bool meta, double kT);
|
||
|
|
|
||
|
|
int main(int argc, char *argv[])
|
||
|
|
{
|
||
|
|
char *data_file;
|
||
|
|
char *out_file;
|
||
|
|
unsigned int step, nsteps, total, out_freq;
|
||
|
|
int *pos, *dpos, *newpos;
|
||
|
|
unsigned int *histogram;
|
||
|
|
const double *grad, *newgrad;
|
||
|
|
unsigned int offset, newoffset;
|
||
|
|
int not_accepted;
|
||
|
|
double dA;
|
||
|
|
double temp;
|
||
|
|
double mbeta;
|
||
|
|
bool meta;
|
||
|
|
double hill, hill_fact, hill_min;
|
||
|
|
double rmsd, rmsd_old, rmsd_rel_change, convergence_limit;
|
||
|
|
bool converged;
|
||
|
|
unsigned int scale_hill_step;
|
||
|
|
|
||
|
|
// Setting default values
|
||
|
|
nsteps = 0;
|
||
|
|
temp = 500;
|
||
|
|
meta = true;
|
||
|
|
hill = 0.01;
|
||
|
|
hill_fact = 0.5;
|
||
|
|
hill_min = 0.0005;
|
||
|
|
|
||
|
|
convergence_limit = -0.001;
|
||
|
|
|
||
|
|
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 << "Version 20110511\n\n";
|
||
|
|
std::cerr << "Syntax: " << argv[0] <<
|
||
|
|
" <filename> [-n <nsteps>] [-t <temp>] [-m [0|1] (metadynamics)]"
|
||
|
|
" [-h <hill_height>] [-f <variable_hill_factor>]\n\n";
|
||
|
|
exit(1);
|
||
|
|
}
|
||
|
|
|
||
|
|
if (meta) {
|
||
|
|
std::cout << "\nUsing metadynamics-style sampling with hill height: " << hill << "\n";
|
||
|
|
if (hill_fact) {
|
||
|
|
std::cout << "Varying hill height by factor " << hill_fact << "\n";
|
||
|
|
}
|
||
|
|
} else {
|
||
|
|
std::cout << "\nUsing unbiased MC sampling\n";
|
||
|
|
}
|
||
|
|
|
||
|
|
if (nsteps) {
|
||
|
|
std::cout << "Sampling " << nsteps << " steps at temperature " << temp << "\n\n";
|
||
|
|
out_freq = nsteps / 10;
|
||
|
|
scale_hill_step = nsteps / 2;
|
||
|
|
converged = true;
|
||
|
|
} else {
|
||
|
|
std::cout << "Sampling until convergence at temperature " << temp << "\n\n";
|
||
|
|
out_freq = 1000000;
|
||
|
|
converged = false;
|
||
|
|
}
|
||
|
|
|
||
|
|
// Inverse temperature in (kcal/mol)-1
|
||
|
|
mbeta = -1 / (0.001987 * temp);
|
||
|
|
|
||
|
|
ABFdata data(data_file);
|
||
|
|
|
||
|
|
if (!nsteps) {
|
||
|
|
scale_hill_step = 2000 * data.scalar_dim;
|
||
|
|
nsteps = 2 * scale_hill_step;
|
||
|
|
std::cout << "Setting minimum number of steps to " << nsteps << "\n";
|
||
|
|
}
|
||
|
|
|
||
|
|
srand(time(NULL));
|
||
|
|
|
||
|
|
pos = new int[data.Nvars];
|
||
|
|
dpos = new int[data.Nvars];
|
||
|
|
newpos = new int[data.Nvars];
|
||
|
|
|
||
|
|
do {
|
||
|
|
for (int i = 0; i < data.Nvars; i++) {
|
||
|
|
pos[i] = rand() % data.sizes[i];
|
||
|
|
}
|
||
|
|
offset = data.offset(pos);
|
||
|
|
} while ( !data.allowed (offset) );
|
||
|
|
|
||
|
|
rmsd = compute_deviation(&data, meta, 0.001987 * temp);
|
||
|
|
std::cout << "\nInitial gradient RMS is " << rmsd << "\n";
|
||
|
|
|
||
|
|
total = 0;
|
||
|
|
for (step = 1; (step <= nsteps || !converged); step++) {
|
||
|
|
|
||
|
|
if ( step % out_freq == 0) {
|
||
|
|
rmsd_old = rmsd;
|
||
|
|
rmsd = compute_deviation(&data, meta, 0.001987 * temp);
|
||
|
|
rmsd_rel_change = (rmsd - rmsd_old) / (rmsd_old * double (out_freq)) * 1000000.0;
|
||
|
|
std::cout << "Step " << step << " ; gradient RMSD is " << rmsd
|
||
|
|
<< " ; relative change per 1M steps " << rmsd_rel_change;
|
||
|
|
if ( rmsd_rel_change > convergence_limit && step >= nsteps ) {
|
||
|
|
converged = true;
|
||
|
|
}
|
||
|
|
|
||
|
|
if (meta && hill_fact && step > scale_hill_step && hill > hill_min ) {
|
||
|
|
hill *= hill_fact;
|
||
|
|
std::cout << " - changing hill height to " << hill << "\n";
|
||
|
|
} else {
|
||
|
|
std::cout << "\n";
|
||
|
|
}
|
||
|
|
}
|
||
|
|
|
||
|
|
offset = data.offset(pos);
|
||
|
|
data.histogram[offset]++;
|
||
|
|
if (meta) {
|
||
|
|
data.bias[offset] += hill;
|
||
|
|
}
|
||
|
|
|
||
|
|
grad = data.gradients + offset * data.Nvars;
|
||
|
|
|
||
|
|
not_accepted = 1;
|
||
|
|
while (not_accepted) {
|
||
|
|
dA = 0.0;
|
||
|
|
total++;
|
||
|
|
for (int i = 0; i < data.Nvars; i++) {
|
||
|
|
dpos[i] = rand() % 3 - 1;
|
||
|
|
newpos[i] = pos[i] + dpos[i];
|
||
|
|
data.wrap(newpos[i], i);
|
||
|
|
if (newpos[i] == pos[i])
|
||
|
|
dpos[i] = 0;
|
||
|
|
|
||
|
|
if (dpos[i]) {
|
||
|
|
dA += grad[i] * dpos[i] * data.widths[i];
|
||
|
|
// usefulness of the interpolation below depends on
|
||
|
|
// where the grid points are for the histogram wrt to the gradients
|
||
|
|
// If done, it has to be done in all directions
|
||
|
|
// the line below is useless
|
||
|
|
//dA += 0.5 * (newgrad[i] + grad[i]) * dpos[i] * data.widths[i];
|
||
|
|
}
|
||
|
|
}
|
||
|
|
|
||
|
|
newoffset = data.offset(newpos);
|
||
|
|
if (meta) {
|
||
|
|
dA += data.bias[newoffset] - data.bias[offset];
|
||
|
|
}
|
||
|
|
|
||
|
|
if ( data.allowed (newoffset) && (((float) rand()) / RAND_MAX < exp(mbeta * dA)) ) {
|
||
|
|
// Accept move
|
||
|
|
for (int i = 0; i < data.Nvars; i++) {
|
||
|
|
pos[i] = newpos[i];
|
||
|
|
not_accepted = 0;
|
||
|
|
}
|
||
|
|
}
|
||
|
|
}
|
||
|
|
}
|
||
|
|
std::cout << "Run " << total << " total iterations; acceptance ratio is "
|
||
|
|
<< double (step) / double (total)
|
||
|
|
<< " ; final gradient RMSD is " << compute_deviation(&data, meta, 0.001987 * temp) << "\n";
|
||
|
|
|
||
|
|
out_file = new char[strlen(data_file) + 8];
|
||
|
|
|
||
|
|
if (meta) {
|
||
|
|
sprintf(out_file, "%s.pmf", data_file);
|
||
|
|
std::cout << "Writing PMF to file " << out_file << "\n";
|
||
|
|
data.write_bias(out_file);
|
||
|
|
}
|
||
|
|
|
||
|
|
// TODO write a PMF for unbiased MC, too...
|
||
|
|
sprintf(out_file, "%s.histo", data_file);
|
||
|
|
std::cout << "Writing sampling histogram to file " << out_file << "\n";
|
||
|
|
data.write_histogram(out_file);
|
||
|
|
|
||
|
|
sprintf(out_file, "%s.est", data_file);
|
||
|
|
std::cout << "Writing estimated FE gradient to file " << out_file << "\n";
|
||
|
|
data.write_field(data.estimate, out_file);
|
||
|
|
|
||
|
|
sprintf(out_file, "%s.dev", data_file);
|
||
|
|
std::cout << "Writing FE gradient deviation to file " << out_file << "\n\n";
|
||
|
|
data.write_field(data.deviation, out_file);
|
||
|
|
|
||
|
|
delete [] pos;
|
||
|
|
delete [] dpos;
|
||
|
|
delete [] newpos;
|
||
|
|
delete [] out_file;
|
||
|
|
exit(0);
|
||
|
|
}
|
||
|
|
|
||
|
|
|
||
|
|
double compute_deviation(ABFdata * data, bool meta, double kT)
|
||
|
|
{
|
||
|
|
// Computing deviation between gradients differentiated from pmf
|
||
|
|
// and input data
|
||
|
|
// NOTE: this is mostly for output, hence NOT performance-critical
|
||
|
|
double *dev = data->deviation;
|
||
|
|
double *est = data->estimate;
|
||
|
|
const double *grad = data->gradients;
|
||
|
|
int *pos, *dpos, *newpos;
|
||
|
|
double rmsd = 0.0;
|
||
|
|
unsigned int offset, newoffset;
|
||
|
|
double sum;
|
||
|
|
int c;
|
||
|
|
bool moved;
|
||
|
|
unsigned int norm = 0; // number of data points summmed
|
||
|
|
|
||
|
|
pos = new int[data->Nvars];
|
||
|
|
dpos = new int[data->Nvars];
|
||
|
|
newpos = new int[data->Nvars];
|
||
|
|
|
||
|
|
for (int i = 0; i < data->Nvars; i++)
|
||
|
|
pos[i] = 0;
|
||
|
|
|
||
|
|
for (offset = 0; offset < data->scalar_dim; offset++) {
|
||
|
|
for (int i = data->Nvars - 1; i > 0; i--) {
|
||
|
|
if (pos[i] == data->sizes[i]) {
|
||
|
|
pos[i] = 0;
|
||
|
|
pos[i - 1]++;
|
||
|
|
}
|
||
|
|
}
|
||
|
|
|
||
|
|
if (data->allowed (offset)) {
|
||
|
|
for (int i = 0; i < data->Nvars; i++)
|
||
|
|
newpos[i] = pos[i];
|
||
|
|
|
||
|
|
for (int i = 0; i < data->Nvars; i++) {
|
||
|
|
est[i] = 0.0;
|
||
|
|
sum = 0.0; // sum of finite differences on two sides (if not on edge of the grid)
|
||
|
|
c = 0; // count of summed values
|
||
|
|
|
||
|
|
newpos[i] = pos[i] - 1;
|
||
|
|
moved = data->wrap(newpos[i], i);
|
||
|
|
newoffset = data->offset(newpos);
|
||
|
|
if ( moved && data->allowed (newoffset) ) {
|
||
|
|
if (meta) {
|
||
|
|
sum = (data->bias[newoffset] - data->bias[offset]) / data->widths[i];
|
||
|
|
c++;
|
||
|
|
} else {
|
||
|
|
if (data->histogram[offset] && data->histogram[newoffset]) {
|
||
|
|
sum = kT * log(double (data->histogram[newoffset]) /
|
||
|
|
double (data->histogram[offset])) / data->widths[i];
|
||
|
|
c++;
|
||
|
|
}
|
||
|
|
}
|
||
|
|
}
|
||
|
|
|
||
|
|
newpos[i] = pos[i] + 1;
|
||
|
|
moved = data->wrap(newpos[i], i);
|
||
|
|
newoffset = data->offset(newpos);
|
||
|
|
if ( moved && data->allowed (newoffset) ) {
|
||
|
|
if (meta) {
|
||
|
|
sum += (data->bias[offset] - data->bias[newoffset]) / data->widths[i];
|
||
|
|
c++;
|
||
|
|
} else {
|
||
|
|
if (data->histogram[offset] && data->histogram[newoffset]) {
|
||
|
|
sum += kT * log(double (data->histogram[offset]) /
|
||
|
|
double (data->histogram[newoffset])) / data->widths[i];
|
||
|
|
c++;
|
||
|
|
}
|
||
|
|
}
|
||
|
|
}
|
||
|
|
|
||
|
|
newpos[i] = pos[i]; // Go back to initial position for next dimension
|
||
|
|
|
||
|
|
est[i] = (c ? sum/double(c) : 0.0);
|
||
|
|
dev[i] = grad[i] - est[i];
|
||
|
|
rmsd += dev[i] * dev[i];
|
||
|
|
norm++;
|
||
|
|
}
|
||
|
|
}
|
||
|
|
|
||
|
|
pos[data->Nvars - 1]++; // move on to next point
|
||
|
|
est += data->Nvars;
|
||
|
|
dev += data->Nvars;
|
||
|
|
grad += data->Nvars;
|
||
|
|
}
|
||
|
|
|
||
|
|
delete [] pos;
|
||
|
|
delete [] newpos;
|
||
|
|
delete [] dpos;
|
||
|
|
|
||
|
|
return sqrt(rmsd / norm);
|
||
|
|
}
|
||
|
|
|
||
|
|
|
||
|
|
char *parse_cl(int argc, char *argv[], unsigned int *nsteps, double *temp,
|
||
|
|
bool * meta, double *hill, double *hill_fact)
|
||
|
|
{
|
||
|
|
char *filename = NULL;
|
||
|
|
float f_temp, f_hill;
|
||
|
|
int meta_int;
|
||
|
|
|
||
|
|
// getting default value for the integer
|
||
|
|
meta_int = (*meta ? 1 : 0);
|
||
|
|
|
||
|
|
// "Syntax: " << argv[0] << " <filename> [-n <nsteps>] [-t <temp>] [-m [0|1] (metadynamics)] [-h <hill_height>]\n";
|
||
|
|
if (argc < 2) {
|
||
|
|
return NULL;
|
||
|
|
}
|
||
|
|
|
||
|
|
for (int i = 2; i + 1 < argc; i += 2) {
|
||
|
|
if (argv[i][0] != '-') {
|
||
|
|
return NULL;
|
||
|
|
}
|
||
|
|
switch (argv[i][1]) {
|
||
|
|
case 'n':
|
||
|
|
if (sscanf(argv[i + 1], "%u", nsteps) != 1)
|
||
|
|
return NULL;
|
||
|
|
break;
|
||
|
|
case 't':
|
||
|
|
if (sscanf(argv[i + 1], "%lf", temp) != 1)
|
||
|
|
return NULL;
|
||
|
|
break;
|
||
|
|
case 'm':
|
||
|
|
if (sscanf(argv[i + 1], "%u", &meta_int) != 1)
|
||
|
|
return NULL;
|
||
|
|
break;
|
||
|
|
case 'h':
|
||
|
|
if (sscanf(argv[i + 1], "%lf", hill) != 1)
|
||
|
|
return NULL;
|
||
|
|
break;
|
||
|
|
case 'f':
|
||
|
|
if (sscanf(argv[i + 1], "%lf", hill_fact) != 1)
|
||
|
|
return NULL;
|
||
|
|
break;
|
||
|
|
default:
|
||
|
|
return NULL;
|
||
|
|
}
|
||
|
|
}
|
||
|
|
|
||
|
|
*meta = (meta_int != 0);
|
||
|
|
return argv[1];
|
||
|
|
}
|